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ABSTRACT  OF  DISSERTATION 


A  DYNAMICAL  THEORY  FOR  HURRICANE  SPIRAL  BANDS 


A  new  theory  is  presented  for  the  formation  of  two  classes  of  hurricane  spiral  bands,  inner 
and  outer.  Inner  bands  are  defined  as  those  which  form  within  500  km  of  the  hurricane 
center  and  owe  their  existence  to  asymmetries  in  the  PV  field.  Outer  bands  are  defined 
as  those  which  form  outside  a  radius  of  500  km  from  the  hurricane  center  and  owe  their 
existence  to  the  breakdown  of  the  intertropical  convergence  zone  (ITCZ).  Inner  bands  are 
discussed  in  the  context  of  Rossby  wave  breaking  and  vortex  merger.  A  shallow  water, 
normal  mode  spectral  model  on  an  /-plane  is  used  to  study  the  banding  process.  Inner 
bands  are  found  to  form  when  the  potential  vorticity  (PV)  field  is  initially  asymmetric  or 
when  a  symmetric  vortex  moves  sufficiently  close  to  another  region  of  high  PV  air.  When 
the  PV  field  is  initially  asymmetric,  wave  breaking  occurs  and  filaments  of  high  PV  are 
ejected  out  to  large  radii  in  the  form  of  spiral  bands.  When  a  symmetric  vortex  moves  close 
to  a  region  high  PV  air,  the  vortex  can  merge  with  the  region  of  PV  and  advect  the  region 
of  PV  towards  the  center  of  the  vortex  in  the  form  of  a  spiral  band.  Because  the  basis 
functions  of  the  model  are  the  normal  modes  of  the  system,  we  can  easily  partition  the  the 
total  flow  into  its  rotational  and  gravitational  components.  The  bands  are  found  to  project 
almost  entirely  onto  the  rotational  modes.  Observations,  however,  show  that  spiral  bands 
are  also  associated  with  strong  convergence,  thus  suggesting  gravitational  modes  play  an 
important  role  in  the  formation  of  spiral  bands.  We  therefore  argue  that  spiral  bands  are 
slow  manifold  phenomena  which  project  onto  both  the  rotational  and  gravitational  modes 
in  such  a  way  that  transient  gravity  waves  are  minimized.  In  addition  to  the  banding 
process,  the  PV  field  was  also  found  to  undergo  a  symmetrization  process.  Langrangian 
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studies  indicate  that  the  higher  PV  air  is  drawn  into  the  center  while  the  lower  PV  air 
appears  to  rotate  about  the  axis,  thus  making  the  PV  field  more  symmetric.  The  wave- 
activity  and  its  flux  are  also  calculated.  These  results  indicate  that  the  bands  are  regions 
of  high  wave-activity  and  therefore  regions  of  high  wave  amplitude.  Experiments  were 
also  performed  in  which  asymmetric  heating  was  added  to  the  model.  When  heating  was 
added,  bands  in  the  PV  field  were  also  formed,  but  in  addition,  the  bands  were  associated 
with  deviations  in  the  geopotential. 

Outer  bands  were  also  studied  with  the  model.  The  model  was  initialized  with  various 
shaped  strips  of  high  PV  air.  These  strips  of  PV  satisfied  the  necessary  condition  for 
barotropic  instability.  When  the  strips  were  perturbed  with  random  noise,  the  PV  began 
to  pool  with  thin  filaments  of  PV  connecting  the  pools.  These  filaments  represent  outer 
spiral  bands  which  are  observed  in  nature.  Results  suggest  that  this  process  can  occur 
anywhere  ITCZ  convection  is  intense. 

The  stability  of  spiral  bands  is  also  considered  using  the  model.  The  bands  are 
idealized  as  annular  regions  of  high  PV  air.  When  a  vortex  of  twice  the  vorticity  of  the 
band  is  introduced,  enough  adverse  shear  is  provided  to  stabilize  the  bands  for  well  over  64 
hours.  When  no  adverse  shear  is  present,  instabilities  in  the  bands  rapidly  grow  and  bands 
are  completely  broken  by  48  horns.  These  experiments  also  suggest  that  polygonal  eyewalls 
and  meso-vortices  may  be  partially  understood  using  barotropic  instability  arguments. 

Thomas  A.  Guinn 

Department  of  Atmospheric  Science 

Colorado  State  University 

Fort  Collins,  Colorado  80523 

Summer  1992 
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Chapter  1 


INTRODUCTION 

One  of  the  most  unmistakable  features  of  tropical  cyclones  is  their  persistent  spiral 
convective  bands.  The  first  known  radar  depiction  of  a  tropical  cyclone  was  taken  at  the 
U.S.  Naval  Air  Station,  Lakehurst,  New  Jersey,  in  September,  1944,  and  the  first  published 
radar  images  of  spiral  bands  were  presented  by  Maynard  (1945).  The  curious  nature  of 
spiral  bands  was  evidenced  in  Maynard’s  paper  when  he  stated, 

“The  echoes  from  this  type  of  weather  phenomena  [tropical  cyclones]  are 
very  easily  identified.  The  eye  of  the  storm  which  appears  as  a  dark  area  is 
surrounded  by  curved  bands  of  echoes  with  feathered  edges  and  trailing  wisps. 

Even  if  the  range  is  not  great  enough  to  encompass  the  eye  of  the  storm,  the 
definite  whorls  of  the  pattern  are  unmistakable.”  Commander  R.H.  Maynard, 

USN,  1945. 

The  curious  nature  of  spiral  bands,  however,  is  not  unique  to  hurricanes.  Astrophysicists 
have  also  struggled  to  explain  the  spiral  nature  of  galaxies,  such  as  that  of  M51.  Numerous 
theories  have  been  proposed  which  attempt  to  explain  this  phenomena  (Toomre,  1977), 
but  a  consensus  has  eluded  the  science.  In  a  similar  fashion,  meteorologists  have  struggled 
to  explain  the  existence  spiral  hurricane  bands,  with  no  true  consensus  being  established 
(although  the  gravity  wave  theory  is  conventionally  accepted).  In  this  study  we  propose 
another  theory  for  hurricane  spiral-  band  formation  in  hopes  that  a  consensus  can  be 
obtained. 

Numerous  meteorological  studies  of  spiral  bands,  both  observational  and  theoretical, 
have  been  accomplished.  Some  of  the  significant  results  of  these  studies,  as  well  as  the 
scope  of  the  present  study,  are  presented  below. 
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1.1  Observed  Features  of  Hurricane  Spiral  Bands 

The  first  published  paper  on  the  observed  structure  of  hurricane  rainbands  was  pre¬ 
sented  by  Wexler  (1947).  He  noted  that  the  passage  of  a  band  was  accompanied  by  heavy 
rain,  a  wind  shift  and  increase  in  wind  speed,  a  large  temperature  fall,  and  a  temporary 
dip  and  recovery  in  the  pressure  of  about  1  mb.  Similar  patterns  were  also  found  by  Ligda 
(1955).  Tatehira  (1961),  in  a  study  of  Typhoon  Helen  (1958),  observed  similar  trends  in 
the  wind  and  mass  fields  but  found  no  significant  change  in  the  temperature  field.  Tatehira 
also  made  calculations  of  the  horizontal  divergence  and  the  vertical  component  of  vorticity. 
These  calculations  showed  the  inner  edge  (closest  the  center)  of  the  band  was  associated 
with  a  convergence  zone  with  a  strong  divergence  zone  near  the  rear  of  the  rainband.  No 
systematic  distribution  of  the  vertical  component  of  vorticity  was  found  however. 

In  more  recent  studies,  Powell  (1990)  and  Ryan  et  ai.  (1992)  have  been  able  to  obtain 
more  extensive  wind  data  through  use  of  airborne  Doppler  radar.  In  both  studies,  wind 
vectors  were  obtained  using  a  pseudo-dual  Doppler  technique.  This  technique  involves 
flying  the  aircraft  in  zigzag  patterns  across  the  bands  to  obtain  radar  scans  in  orthogonal 
directions.  The  winds  were  assumed  to  be  steady  for  the  time  required  to  complete  two 
orthogonad  scans,  thus  allowing  two-dimensional  wind  vectors  to  be  calculated.  The  use 
of  this  technique  has  allowed  for  the  calculation  of  reasonably  accurate  and  dense  wind 
fields.  The  most  pertinent  results  to  the  present  study  sure  the  calculations  of  the  relative 
vorticity  and  horizontal  divergence. 

Powell  (1990)  presented  hurricane  boundary  layer  results  from  two  tropical  cyclones, 
Josephine  (1984)  and  Earl  (1986).  In  each  of  these  cases,  the  bands  were  stationary  with 
respect  to  the  storm.  At  the  time  of  Josephine’s  analysis,  the  maximum  flight  level  winds 
were  50  ms-1  and  the  minimum  sea-level  pressure  was  970  mb.  At  the  time  of  Earl’s 
analysis,  the  maximum  flight  level  Winds  were  also  50  ms-1  and  the  pressure  was  983 
mb.  Cross-band  fields  from  a  relatively  narrow  band  (~20  km)  of  Josephine  showed  a 
region  of  high  relative  vorticity  associated  with  the  band  with  a  maximum  of  2.5  xlO-3 
s'1  at  1500  m.  In  addition,  the  inner  edge  of  the  band  was  associated  with  a  maximum 
in  convergence  of  2.0xl0-3  s~l  at  500  m.  The  band  examined  in  Hurricane  Earl,  which 
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was  slightly  wider  (~30  km),  exhibited  the  same  general  characteristics  as  Josephine  but 
the  maximums  were  slightly  lower  in  magnitude.  The  high  relative  vorticity  associated 
with  the  bands  in  both  cases  was  due  to  strong  cyclonic  shear  across  the  band  with  a 
maximum  of  the  along-band  wind  component  near  the  outer  edge  of  the  band.  Other 
dynamical  features  which  were  common  to  both  Josephine  and  Earl  were  upward  vertical 
velocity  and  convergence  maxima  on  the  inner  side  of  the  band  axis  and  D-value  minima 
found  at  the  band  axis.  The  D-value  is  the  difference  between  the  radar  altitude  and 
the  pressure  altitude  with  reference  to  the  standard  atmosphere.  It  is  therefore  analogous 
to  the  deviation  geopotential.  In  Fig.  1.1,  cross-band  fields  of  Hurricane  Josephine  from 
Powell  (1990)  are  presented.  This  figure  shows  vertical  cross  sections  based  on  two  passes 
through  the  band,  one  at  500  m  and  one  at  1500  m.  The  dashed  box  in  panel  (a)  represents 
a  cross  section  in  the  x-z  plane  (i.e.  coming  out  of  the  page). 

Ryan  et  aJ.  (1992)  examined  a  relatively  wide  band  associated  with  Tropical  Cyclone 
Irma  (1987)  which  formed  over  the  Gulf  of  Carpenteria.  The  band  examined  in  their  study 
was  relatively  wide  (~  80  km).  At  the  time  of  the  analysis,  the  observations  indicated  Irma 
was  at  least  a  strong  depression.  It  was  declared  to  be  a  named  tropical  storm  roughly 
15  hours  later,  achieving  a  minimum  pressure  of  978  mb.  Although  Ryan  et  ad.  did  not 
explicitly  calculate  the  divergence  and  vorticity,  they  did  calculate  both  the  along-band 
and  cross-band  components  of  the  the  wind  fields.  These  calculations  are  shown  in  Fig. 
1.2  along  with  a  radar  composite  of  the  storm.  To  estimate  the  vorticity  and  divergence 
fields,  we  used  the  same  calculation  as  Powell  (1990),  i.e. 


div  ~ 


dVcro 

dx 


(1.1) 


and 


c« 


dValo 
dx  ' 


(1.2) 


where  x  is  the  cross-band  distance,  and  Vcro  and  V^o  are  the  cross-band  and  along-band 
components  of  the  wind,  respectively.  This  allowed  us  to  conclude  that  the  band  had 
similar  vorticity  and  divergence  patterns  to  those  of  Josephine  and  Earl.  They  were, 
however,  roughly  an  order  of  magnitude  weaker.  This  can  be  attributed  to  the  fact  that 
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Figure  1.1:  Cross-band  dynamical  fields  of  Hurricane  Josephine  (1984).  (a)  40x40  km 
lower  fuselage  radar  composite,  (b)  band-relative  windspeed,  (c)  along-band  wind  com¬ 
ponent,  (d)  cross-band  wind  component,  (e)  divergence,  (f)  vertical  wind,  (g)  relative 
vorticity  and  (h)  D- value  perturbation.  (From  Powell,  1990.) 
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Figure  1.2:  Cross- band  dynamical  fields  from  Tropical  Depression  Irma  (1987).  (a)  Re¬ 
flectivity  pattern  from  Gove  and  P-3  radars  showing  the  major  rainband  and  its  relation 
to  the  center  of  the  incipient  cyclone.  The  heavy  line  locates  the  aircraft  tracks  used  for 
cross  sections,  (b)  Composite  analysis  of  the  cross-band  wind  component  (ms-1).  The 
rainband  center  is  near  70  km  and  the  dashed  lines  show  the  approximate  location  of  the 
rainband  edges,  (c)  The  same  as  (b)  but  for  the  along-band  component.  (From  Ryan  et 
a/.,  1992.) 
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Figure  1.3:  Wind  (knots)  and  D-value  (geopotential  meters)  for  the  550  mb  surface.  The 
trough  in  D-value  is  clearly  associated  with  the  rainband.  (From  Ryan  et  a 1.,  1992.) 

Irma  was  not  yet  a  tropical  cyclone  at  the  time  of  the  analysis  and  that  Iona  was  a  much 
wider  band.  Even  if  Irma  had  achieved  the  cross-band  wind  shears  of  Josephine  and 
Earl,  the  divergence  and  vorticity  would  still  be  much  weaker  because  of  the  width  over 
which  the  shear  occurred.  Ryan  et  a J.  also  demonstrated  that  the  bands  were  associated 
with  minima  in  the  D-value.  Figure  1.3  shows  the  D-value  field  and  winds  at  550  mb 
associated  with  TVopical  Depression  Irma.  Similar  to  Powell’s  results,  comparisons  with 
Fig.  1.2a  show  that  the  trough  in  D-value  also  appears  to  be  nearly  centered  about  the 
band  axis  in  Irma  as  well.  The  significance  of  all  these  bands  displaying  relative  maxima 
in  both  the  divergence  and  vorticity  fields  will  be  made  clear  shortly. 

1.2  Theories  of  Spiral  Band  Formation 

Since  the  discovery  of  spiral  rainbands,  many  theories  have  been  presented  to  explain 
their  existence.  Fletcher  (1945)  first  theorized  that  the  spiral  bands  from  tropical  cyclones 
which  had  formed  on  the  intertropical  convergence  zone  (1TCZ)  were  produced  by  “indi¬ 
vidual  convergence  lines”  which  had  been  coiled  in  the  the  centers  of  the  storms.  Wexler 
(1947)  hypothesized  that  the  phenomena  wasn’t  unique  to  the  ITCZ,  but  rather  bands 
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could  occur  any  time  a  vortex  came  in  contact  with  existing  “cloud  streets.”  These  cloud 
streets  would  then  become  curved  and  the  convergence  associated  with  the  vortex  would 
cause  the  bands  to  spiral  inward  towards  the  center.  The  advent  of  satellite  meteorology, 
however,  has  negated  these  theories. 

Tepper  (1958)  later  related  hurricanes  to  mid-latitude  squall  lines.  He  proposed  the 
bands  were  gravity  waves  which  propagated  out  from  the  center.  If  there  existed  a  stable 
layer  aloft,  and  “something”  impressed  an  inflow  on  the  circulation  of  the  vortex,  gravity 
waves  would  be  induced  near  the  center  of  the  vortex  and  propagate  outwards  at  the 
interface  of  the  stable  layer.  Spiral  regions  of  convection  could  then  be  generated  if  the 
height  of  the  stable  layer  varied  azimuthally  or  if  the  tropical  cyclone  was  observed  to 
move  rapidly.  Abdullah  (1966)  also  related  hurricane  bands  to  gravity  waves  at  the  upper 
surface  of  the  lower  layer.  He  suggested  that  gravity  waves  were  generated  by  a  fresh 
surge  of  air  at  the  periphery  of  the  storm  that  propagated  towards  the  storm  center. 
Upon  reaching  the  eye,  the  waves  would  break  causing  turbulence  and  the  release  of 
latent  instability.  He  argued  that  this  would  manifest  itself  in  the  form  of  a  squall  which 
could  then  grow  and  propagate  back  outwards,  taking  on  a  spiral  shape.  He  presented  a 
simple  two  layer  model  to  justify  his  theory.  Kurihara  (1976)  suggested  the  bands  were 
outward  propagating  gravity  waves  which  were  modified  slightly  by  rotation  (i.e.  gravity- 
inertia  waves).  He  also  noted  that  bands  could  also  be  produced  in  his  six  layer  model  by 
two  other  types  of  waves  which  propagated  inwards,  one  of  which  was  simply  the  inward 
propagating  gravity-inertia  wave.  The  third  type  of  wave,  however,  was  found  to  have 
features  of  a  geostrophic  mode.  He  noted,  however,  that  the  third  type  had  not  been 
observed  in  nature.  In  contrast,  Willoughby  (1978)  suggested  the  inward  propagating 
gravity-inertia  waves  were  most  likely  to  occur  in  nature  and  therefore  most  likely  to  be 
responsible  for  producing  the  spiral  bands.  If  spiral  bands  in  nature  were  due  solely  to 
gravity-inertia  modes,  we  could  not  explain  the  observed  relative  vorticity  maxima  which 
are  associated  with  observed  spiral  bands  (Powell,  1990;  Ryan,  1992).  It  therefore  seems 
likely  that  these  bands  are  not  strictly  gravity  waves. 

Although  the  gravity  wave  theory  is  presently  accepted  by  convention,  other  theories 
have  also  been  presented  which  attempt  to  explain  the  formation  of  spiral  bands.  Faller 
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(1961)  suggested  that  spiral  bands  were  formed  as  a  result  of  boundary  layer  instabilities. 
His  argument  was  based  on  laboratory  experiments  of  a  fluid  in  a  rotating  tank.  When  the 
Reynold’s  number  of  the  flow  exceeded  145,  roll  vortices  appeared.  These  roll  vortices  were 
observed  to  have  crossing  angles  similar  to  those  of  hurricane  spiral  bands.  Fung  (1977) 
also  considered  boundary  layer  instability  as  a  possible  mechanism  for  hurricane  band 
formation.  In  her  study  she  considered  Rayleigh  instability  as  the  mechanism  responsible 
for  band  formation.  She  noted  that  Rayleigh  instability  occurred  in  flows  in  which  the 
Reynold’s  number  exceeded  110.  The  roll  vortices  formed  by  Rayleigh  instability  were 
calculated  to  be  oriented  14°  to  the  left  of  the  geostrophic  wind.  The  existence  of  these 
roll  vortices  also  depended  highly  on  the  existence  of  an  inflection  point  in  the  radial 
wind  profile,  which  is  easily  satisfied  since  the  radial  wind  is  observed  to  have  a  maximum 
in  the  boundary  layer.  The  results  of  Fung’s  stability  analysis  indicated  that  the  most 
unstable  mode  had  four  radial  arms  and  radial  wavelengths  of  about  20  km  near  the 
eyewall  and  about  50  km  at  a  radius  of  300  km.  The  spiral  pattern  was  observed  to  be 
nearly  stationary.  In  the  present  study,  we  obtain  spiral  bands  without  the  presence  of  a 
boundary  layer.  Although  there  is  no  direct  conflict  between  the  theory  presented  here 
and  boundary  layer  instability  theories,  the  theory  presented  in  this  study  helps  explain 
spiral  bands  in  a  much  simpler  fashion. 

The  final  and  most  overlooked  theory  of  spiral  bands  is  that  of  Rossby  waves.  Mac¬ 
Donald  (1968)  argued  that  tha  spiral  bands  in  hurricanes  bore  a  resemblance  to  troughs 
in  the  general  circulation  of  the  atmosphere.  The  clouds  associated  with  these  troughs 
would  appear  to  spiral  outwards  when  seen  from  an  extra-terrestrial  viewpoint.  An  artists 
conception  of  how  this  might  appear  is  shown  in  Fig.  1.4.  MacDonald’s  analogy  received 
little  attention  because  at  the  time  of  its  publication  (and  well  after)  there  were  no  obser¬ 
vational  studies  which  showed  bands  to  be  associated  with  regions  of  high  relative  vorticity 
as  would  be  expected  if  they  were  Rossby  waves  (e.g.,  Anthes,  1982).  In  light  of  recent 
observations  (e.g.,  Powell,  1990;  Ryan  et  a J.,  1992),  MacDonald’s  Rossby  wave  theory 
deserves  more  attention. 
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Figure  1.4:  An  artist’s  conception  of  how  clouds  associated  with  tilted  troughs- in  the 
troposphere  might  appear  from  an  extra  terrestrial  viewpoint.  (From  MacDonald,  1968.) 

1.3  The  Scope  of  the  Present  Study 

In  this  study,  we  propose  that  spiral  bands  are  neither  entirely  gravity  wave  phenom¬ 
ena  nor  entirely  Rossby  wave  phenomena.  Observationally,  this  seems  likely  since  bands 
are  associated-  with  both  horizontal  convergehce/divergence  and  relative  vorticity.  We 
therefore  propose  that  spiral  bands  are  slow  manifold  phenomena.  That  is,  they  project 
onto  both  the  gravity-inertia  modes  as  well  as  the  Rossby  (rotational)  modes  in  such  a 
way  that  transient  gravity  waves  are  minimized.  The  formation  of  these  spiral  bands  can 
easily  be  understood  through  use  of  wave  breaking  and  vortex  merger  arguments.  Previ¬ 
ously,  wave  breaking  arguments  have  only  been  used  to  study  stratospheric  dynamics,  e.g., 
McIntyre  and  Palmer  (1984),  Juckes  and  McIntyre  (1987),  and  Polvani  and  Plumb  (1992). 
Several  wave  breaking  studies  (e.g.,  Dritschel,  1989,1992;  Polvani  and  Plumb,  1992)  have 
been  accomplished  using  contour  dynamical  methods.  In  contour  dynamics,  the  positions 
of  fluid  particles  which  define  a  material  contour  are  predicted.  This  method  provides 
an  elegant  means  by  which  the  morphology  of  contour  shapes  cam  be  examined.  These 
methods,  however,  require  the  absence  of  diabatic  or  frictional  effects  since  these  destroy 
the  ability  to  define  a  materially  conserved  variable.  In  the  present  model,  experiments 
are  performed  using  a  shallow  water  normal  mode,  spectral  model  on  am  /-plane.  Since 
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the  basis  functions  of  a  normal  mode  model  are  the  solutions  to  the  linear  system,  they 
provide  a  natural  means  with  which  to  partition  the  gravity-inertia  modes  from  the  rota¬ 
tional  modes.  Because  of  this  feature,  we  can  see  directly  the  part  of  the  solution  which 
is  due  to  the  gravity-inertia  modes  and  the  part  which  is  due  to  the  rotational  modes. 

We  provide  dynamical  explanations  for  two  types  of  hurricane  bands,  inner  and  outer. 
Take  note,  however,  that  our  definitions  of  inner  and  outer  bands  are  not  those  which 
are  conventionally  accepted.  We  define  inner  bands  as  those  bands  which  typically  form 
within  500  km  of  the  vortex  center  and  owe  their  existence  to  asymmetries  in  the  potential 
vorticity  field.  Outer  bands,  on  the  other  hand,  are  those  which  typically  form  outside  of 
500  km  from  the  center  and  owe  their  existence  to  the  breakdown  of  the  ITCZ.  Both  wave 
breaking  and  vortex  merger  processes  are  suggested  to  be  responsible  for  the  formation  of 
inner  bands,  while  barotropic  instability  is  suggested  to  be  primarily  responsible  for  outer 
bands. 

The  outline  of  the  study  is  as  follows.  In  Chapter  2,  we  present  the  numerical  model 
which  is  used  for  all  simulations.  We  include  in  this  chapter  the  normal  mode  transform 
which  allows  us  to  partition  the  the  gravity-inertia  modes  horn  the  geostrophic  (rotational) 
modes.  In  Chapter  3,  we  examine  the  formation  of  inner  bands  by  both  wave  breaking 
and  vortex  merger.  Wave-activity  and  flux  of  wave- activity  are  calculated  to  help  with 
diagnostics.  Another  diagnostic  tool  used  in  this  chapter  is  passive  tracer  experiments. 
This  tool  allows  us  to  examine  Lagrangian  trajectories  of  the  fluid.  Also  in  Chapter  3, 
we  consider  the  evolution  of  a  vortex  which  is  heated  asymmetrically.  In  Chapter  4, 
we  consider  the  formation  of  outer  bands  by  the  breakdown  of  the  ITCZ.  Several  general 
stability  theorems  are  presented  to  aid  in  understanding  the  instabilities  involved.  Normal 
mode  stability  analysis  is  also  used  to  gain  insight  into  the  unstable  patterns  and  growth 
rates.  Five  numerical  simulations  are  presented  which  demonstrate  the  breakdown  of  the 
ITCZ  and  the  formation  of  outer  bands.  In  Chapter  5,  we  consider  the  stability  of  spiral 
bands.  Again  normal  mode  stability  analysis  is  used  to  predict  unstable  growth  rates  and 
patterns.  These  experiments  also  lead  us  to  a  proposed  theory  for  polygonal  eyewalls  and 
meso- vortices.  Finally,  in  Chapter  6,  we  present  a  brief  summary  and  some  additional 
conclusions. 
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It  needs  to  be  noted  that  the  terms  “hurricane”  and  “tropical  cyclone”  are  used 
interchangeably  throughout  the  text.  “Hurricane”  never  refers  specifically  to  those  tropical 
cyclones  which  form  east  of  the  dateline,  unless  we  refer  to  a  specific  named  storm. 


Chapter  2 


THE  NORMAL  MODE  SPECTRAL  MODEL 

In  this  chapter,  we  present  the  normal  mode  shallow  water  spectral  model  which 
is  used  in  all  the  numerical  simulations.  Included  are  discussions  of  the  normal  mode 
transform,  the  calculation  of  the  nonlinear  terms  by  the  transform  method,  the  time 
discretization,  and  the  calculation  of  derivatives.  In  the  final  section  we  discuss  some  of 
the  advantages  and  disadvantages  of  normal  mode  spectral  models  as  compared  with  finite 
difference  models  and  other  spectral  models. 

2.1  The  Governing  Equations 

Using  Cartestian  coordinates,  the  shallow  water  equations  on  an  /-plane  can  be  writ¬ 
ten  in  rotational  form  as 

tr  - vU + ° + b  [* + i(“2 + ”2)1 = F 

+  “(/  +  0  +  ^  [*  +  £(“2  +  »2)]  =  <j 

d<f>  ,  2 /Su  ,  dv\  ,  d(u<f>)  d(v<f>) 

at  + *  \ai  +  a} )  +  ~dT  +  ST  = 

where  u  and  v  are,  respectively,  the  east-west  and  north-south  components  of  the  wind 
field.  In  should  be  noted,  however,  that  their  orientation  is  completely  arbitrary  on  an 
/-plane.  The  geopotential,  is  defined  as  gh  with  g  and  h  being,  respectively,  gravity  and 
the  deviation  of  the  fluid  depth  from  its  mean  value,  H.  The  variable,  c,  is  defined  as  >/gH. 
Frictional  affects  are  represented  by  F  and  G  while  diabatic  affects  are  represented  by  S. 
Lastly,  the  variable  £  represents  the  relative  vorticity  which  is  given  as  (dv/dx  -  du/dy). 


(2.1) 

(2.2) 

(2.3) 
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A  more  convenient  way  to  write  (2.1)-(2.3)  is  in  vector  form.  This  is  done  by  defining 
W  =  [«,t >,<f>/c\T  and  rewriting  (2.1)-(2.3)  as 


^-  +  £W  =  F  +  N, 

where  the  linear  operator  £  is  defined  as 

f  0  “/  > 

^  =  /  0 

<  ch  cvi  0  > 

nonlinear  and  forcing  terms  are  written  as 


(2-4) 


(2.5) 


and 


F  =  (F,  G,  S/c)T , 


(2.7) 


respectively.  With  the  governing  equations  written  in  this  form,  we  can  now  begin  the 
normal  mode  transform. 


2.2  The  Normal  Mode  Transform 


We  begin  the  normal  mode  transform  by  first  defining  an  appropriate  inner  product. 
If  we  consider  our  model  domain  to  be  periodic  in  both  the  x  and  y  directions  with  periods 
Lx  and  Ly  respectively,  we  can  define  our  complex  inner  product  to  be 

W,0=  tV  fL'  r^oCo+^lVx+^dxdy  (2.8) 

JO  JO 

where  V>  and  (  are  three-component  complex  vectors  which  are  periodic  in  x  and  y  with 
periods  Lx  and  Ly  respectively  and  the  asterisk  denotes  the  complex  conjugate.  As  dis¬ 
cussed  in  Appendix  B,  £  is  skew-Hermeti&n  with  respect  to  the  inner  product  (2.8).  By 
noting  that  the  eigenvalues  of  skew-Hermetian  operators  are  pure  imaginary,  we  can  write 

LKlclq  =  IVklqKklq,  (2.9) 

where  i/y,  are  the  eigenvalues  of  the  matrix  operator  £  and  K klq  are  the  corresponding 
vector  eigenfunctions  or  transformation  kernels.  The  subscripts  k,  l,  and  q  are  presently 
undefined.  They  will  become  clear  once  (2.9)  is  solved. 
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Equation  (2.9)  can  easily  be  solved  in  the  usual  manner  by  first  assuming  that  the 
three  components  of  are  proportional  to  e^kx+ly\  where  i  =  y/—l.  By  assuming 
this,  all  x  derivatives  are  replaced  by  tk  and  all  y  derivatives  are  replaced  by  il.  The 
matrix  operator  C  is  then  reduced  to  a  simple  coefficient  matrix.  In  order  for  a  non-trivial 
solution  to  exist,  the  determinant  of  this  coefficient  matrix  must  vanish.  We  find  that  the 
determinant  will  vanish  provided  the  following  dispersion  relation  holds, 


VklqV'hn  ~  ib2  +  f2)]  =  0, 


(2.10) 


where 


fi2  =  ^(k2  +  l2). 


(2.11) 


Thus,  our  eigenvalues  are 


0  q  =  0  (geostrophic  modes) 

i 'kiq  =  <  —  (/2  +  /x2)a  9=1  (gravity  -  inertia  mode) 

(/2  +  /i2)  j  9  =  2  (gravity  —  inertia  mode) 


(2.12) 


With  the  use  of  (2.12),  the  transformation  kernels,  which  are  simply  the  eigenvectors  of 
C,  can  easily  be  determined.  They  are 


(2.13) 


(  ~id  ^ 

£  I  ick  e»(fc*+<v)  9  =  0 

A  /  /  v 

/  c (if  l  -  vk)  \ 

K*,,=  <  Tfc  -cH  +  i/k)  9  =  1 

/  c{vk  +  ifl)  \ 

^  c(vl-ifk)  9  =  2 


where  i *  =  [p2  +  /*]*. 

The  definitions  of  the  subscripts  should  now  be  clear.  The  subscripts  k  and  l  represent 
wave  numbers  in  the  x  and  y  directions  respectively,  while  the  subscript  q  refers  to  a 
particular  root  of  (2.10).  Physically,  q  distinguishes  between  the  fast  and  slow  modes. 
That  is,  9  =  0  denotes  the  geostrophic  (slow)  modes  while  9  =  1  or  2  denotes  the  east  and 


west  propagating  gravity-inertia  (fast)  modes,  respectively. 
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As  discussed  in  Appendix  C,  the  transformation  kernels  in  (2.13)  are  orthonormal 
with  respect  to  the  complex  inner  product  (2.8),  that  is 


(Kklq,Kkniq>) 


k  =  k\  V  =  /  and  <f  =  q 
otherwise 


}• 


Using  this  orthonormality,  we  can  define  our  transform  pair  as 

w(x,» ,t)=  f;  f;  x;w«f(t)K«,(*,tf) 

/=— oo  fc=— oo  9=0 

WrWo(t)  =  (W(x,y,t),Kw,(x,y)), 


(2.14) 


(2.15) 

(2.16) 


where  Wkiq(t)  is  a  scalar,  complex  function  of  time  having  units  of  ms-1.  Equation  (2.15) 
can  then  be  thought  of  as  the  transformation  which  takes  the  dependent  variables  from 
physical  space  to  spectral  space,  and  (2.16)  is  the  inverse  transformation  back  to  physical 
space. 

Finally,  we  use  our  spectral  transforms  (2.15)  and  (2.16)  to  reduce  our  original  system 
of  equations  to  one  set  of  ordinary  differential  equations  for  each  mode.  We  accomplished 
this  by  taking  the  inner  product  (as  defined  by  Eq.  2.8)  of  (2.4)  with  the  transformation 
kernel  K kiq(x,  y)  to  obtain 


J~(W,  K«,)  +  (£W,K«,)  =  (F,K«,)  +  (N,KW*). 


(2.17) 


Because  £  is  skew-Hermetian  with  respect  to  our  inner  product,  we  can  shift  £  off  of  W 
and  onto  the  transformation  kernel  (see  Appendix  B).  Using  (2.9),  we  can  then  eliminate 
£  from  the  problem  which  allows  us  to  write  (2.17)  in  final  form  as 

+  ivkl9Wklq  =  Fklq  +  Nklq.  (2.18) 


Equation  (2.18)  is  the  normal  mode  oscillation  equation.  In  the  absence  of  external  forcings 
or  nonlinear  terms,  each  mode  would  simply  oscillate  at  its  own  frequency.  It  is  interesting 
to  note  that  all  normal  mode  models  can  be  reduced  to  this  one  equation  regardless  of 
geometry.  In  Appendix  A,  we  demonstrate  this  by  deriving  the  normal  mode  oscillation 
equation  in  cylindrical  coordinates. 
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2.3  Computational  Procedures 

In  this  section,  we  discuss  the  computational  procedures  involved  in  numerically  inte¬ 
grating  the  above  system  of  equations.  These  procedures  include  the  numerical  calculation 
of  the  finite  Fourier  and  normal  mode  transforms.  We  also  discuss  the  calculation  of  deriva¬ 
tives  as  well  as  nonlinear  terms.  Lastly,  we  discuss  the  time  differencing  scheme  which 
was  used. 

2.3.1  Computational  Aspects  of  the  Horizontal  Transforms 

In  the  preceding  section,  the  normal  mode  transform  was  described  as  a  single  op¬ 
eration.  While  this  is  a  convenient  way  to  express  the  transform  on  paper,  it  is  not 
convenient  for  computational  purposes.  Computationally,  it  is  much  easier  to  perform  the 
same  normal  mode  transform  as  two  separate  operations.  In  addition,  we  must  truncate 
our  transforms  in  order  to  compute  them  numerically. 

The  first  operation  consists  of  transforming  the  wind  and  geopotential  fields  to  Fourier 
spectral  space.  We  define  our  truncated  Fourier  transform  pair  as 


4>u  =  JJ  H  ]C ^xbVj) exP  |-2,r*  (t  +  J )] 

(2.19) 

V>(*;, yj)  =£  £  V-fc/exp  [2**  (7  +  j]\ . 
l--Lk=-K  L  \  /  J 

(2.20) 

where  the  continuous  variables  have  been  discetized  as 

k;  =  7— fc,  for  k  =  -K,...,K  (2.21) 

l:  —  7 —l,  for  l  =  —L, ...,  L  (2.22) 

Ly 

i\  =  ^j-i  for  i  =  0, ...,  7-1  (2.23) 

Vj  =  j  for  3  =  0, ...,  7-1  (2.24) 

with  7  and  7  representing  the  maximum  number  of  grid  points  on  the  transform  grid  in  the 
x  and  y  directions  respectively,  and  K  and  L  represent  the  maximum  number  of  discrete 
wave  numbers  in  the  x  and  y  directions  repsectively.  Equation  (2.19)  is  the  transform 
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from  physical  to  Fourier  spectral  space  and  (2.20)  is  its  inverse.  In  order  to  speed  the 
model  computations,  fast  Fourier  transforms  (FFTs)  are  used  to  evaluate  (2.19-2.20)  (see 
Appendix  D).  If  the  number  of  gridpoints  used  in  both  model  directions,  I  and  J,  are  both 
powers  of  two,  FFTs  reduce  the  number  of  operations  needed  to  evaluate  either  (2.19)  or 
(2.20)  from  order  72  x  J2  to  order  I  log2 1  x  Jlog2  J  (Conte  and  de  Boor,  1980).  As  an 
example,  if  128  x  128  gridpoints  are  used  to  represent  the  physical  fields,  FFTs  reduce  the 
number  of  operations  needed  to  transform  these  fields  to  or  from  Fourier  spectral  space 
from  order  108  to  order  106.  This  is  an  improvement  of  roughly  two  orders  of  magnitude. 

Once  the  fields  have  been  transformed  to  Fourier  spectral  space,  the  second  operation, 
which  can  be  thought  of  as  the  true  normal  mode  transform,  reduces  to  a  simple  vector 
multiplication  of  the  Fourier  spectral  coefficients  with  the  complex  conjugate  of  the  x  and 
y  independent  part  of  the  transformation  kernel  (2.13)  for  each  mode  type  q. 

The  inverse  transformation  is  the  summation  of  the  spectral  coefficients,  Wkiq  with 
the  x  and  y  independent  parts  of  the  transformation  kernel  for  each  q.  Once  this  is  done, 
the  inverse  Fourier  transform  is  again  performed  via  FFTs. 

2.3.2  Calculation  of  Derivatives 

An  advantage  of  working  with  a  doubly  periodic  domain  is  that  the  transformation 
kernels  are  eigenfunctions  of  the  linear  operators;  that  is,  derivatives  of  the  transformation 
kernels  are  simply  scalar  multiples  of  the  transformation  kernels  themselves.  We  can  use 
this  property  to  show  that  the  spectral  coefficients  of  derivatives  of  our  fields  are  simply 
scalar  multiples  of  the  spectral  coefficients  of  our  fields  themselves.  To  show  this,  consider 
the  one  dimensional  finite  Fourier  transform  of  the  periodic  function,  drp(x)/dx.  We  write 


this  as 


1  [L*  fop(x)  ik 


-1/ 

Lz  Jo 


(2.25) 


If  we  integrate  (2.25)  by  parts,  we  can  shift  the  derivative  operator  off  t/>(x)  to  obtain 


j-  ikj,(x)e-ikxdx  +  [*(*)«-**]*'  . 


(2.26) 
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Because  ip  is  periodic  in  x  the  rightmost  term  in  (2.26)  is  identically  zero.  This  implies 
that 

i  fL‘  =  .1 


cdx  =  ikrfrk, 


(2.27) 


-l 

Lx  Jo  dx 

where  xpk  is  the  kth  Fourier  spectral  coefficient  of  i/>(x).  Although  this  proof  is  only  in  one 
dimension,  it  can  easily  be  generalized  to  two  dimensional  transforms  and  therefore  our 
inner  product  (2.8).  We  can  therefore  write 


=  ^(W(x,y,t),K  «,(*,y)).  (2.28) 

This  is  not  only  true  for  x  derivatives  but  similar  relations  also  hold  true  for  all  linear 
operators,  e.g.  second  order  derivatives,  Laplacians,  etc.  This  is  useful  because  it  requires 
us  to  store  only  one  transformation  kernel.  All  necessary  derivatives  can  then  be  obtained 
through  scalar  multiplications  of  this  single  kernel  given  in  (2.13). 


2.3.3  Evaluation  of  Nonlinear  Terms  Using  the  Transform  Method 


Evaluation  of  the  nonlinear  terms  of  the  model  is  the  most  expensive  computation 
of  the  model.  The  reason  being  that  nonlinear  products  cannot  be  computed  in  spectral 
space  simply  by  multiplying  the  appropriate  Fourier  coefficients.  In  this  model  we  use 
the  transform  method  as  described  by  Orzag  (1970)  and  Eliasen  et  aJ.  (1970)  to  calculate 
the  nonlinear  terms.  The  transform  method  predicts  the  linear  terms  and  calculates 
their  derivatives  in  spectral  space  while  the  nonlinear  terms  are  computed  in  physical 
space.  To  demonstrate  this,  consider  the  normal  mode  oscillation  equation  (2.18).  The 
nonlinear  terms,  Nkiq,  are  necessary  to  predict  the  dependent  variables,  Wkiq,  however 
the  calculation  of  Nktq  depends  on  Wkiq  itself.  To  calculate  Nkiq  at  a  given  timestep,  we 
perform  the  inverse  normal  transform  of  Wktq  at  the  corresponding  time  step  to  obtain 
the  Fourier  coefficients.  While  in  Fourier  spectral  space,  all  appropriate  derivatives  are 
calculated.  We  next  transform  the  fields  and  their  derivatives  to  physical  space  and 
compute  the  nonlinear  products.  The  nonlinear  products  are  then  transformed  back  to 
to  normal  mode  spectral  space  to  obtain  Nkiq •  The  forcing  terms,  Fkiq  are  calculated  in 
a  similar  fashion.  Using  Nkiq  and  Fklq,  a  new  Wkiq  can  be  predicted  from  (2.18)  and  the 
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whole  process  repeats  itself.  It  should  now  be  easy  to  see  why  the  nonlinear  and  forcing 
terms  use  such  a  large  portion  of  the  total  computing  time.  For  each  time  step,  a  complete 
transform  to  and  from  spectral  space  is  required. 

It  is  important  to  note  that  because  our  basis  functions  are  periodic,  the  Fourier 
transforms  can  be  computed  exactly  if  enough  grid  points  are  used  to  discretize  the  fields 
in  physical  space.  This  holds  true  even  for  quadratically  nonlinear  terms.  To  demonstrate 
this,  recall  Krylov  (1962)  where  it  is  shown  that  any  periodic  function  of  period  Lx  with 
degree  no  greater  than  K  can  be  integrated  over  Lx  exactly  by  trapezoidal  quadrature  if  at 
least  K  + 1  grid  points  are  used  to  discretize  the  function.  Now  consider  any  quadratically 
nonlinear  term.  Each  function  which  comprises  the  nonlinear  term  is  of  degree  of  less  than 
or  equal  to  K  and  L;  that  is,  each  function  can  be  Fourier  decomposed  to  include  only 
terms  of  the  form  e'^kx+ly^  where  k  <  K  and  l  <  L.  Since  the  integral  transform  itself  (2.8) 
also  includes  the  term  e‘(fcx+*v)  for  k  <  K  and  l  <  L,  the  highest  degree  of  the  resulting 
function  must  be  3 K  and  3 L,  i.e.  the  resulting  function  takes  the  form  e^3kx+3ly\  This 
means  we  can  evaluate  the  integral  in  (2.8)  exactly  for  up  to  and  including  quadratically 
nonlinear  terms  if  we  use  at  least  3  K  +  1  and  3L  +  1  points  in  the  x  and  y  directions 
respectively.  As  an  example,  if  we  choose  a  128  x  128  point  grid,  we  can  resolve  roughly 
42  waves. 

If  all  terms  in  the  model  were  strictly  linear  or  quadratically  nonlinear,  all  transfor¬ 
mations  could  be  done  exactly.  As  a  result,  no  aliasing  error  would  occur;  that  is,  no 
information  from  higher  order  Fourier  modes  would  be  falsely  projected  onto  the  lower 
order  Fourier  modes.  If  this  were  true,  non-linear  computational  instability  arising  from 
this  false  projection  as  described  by  Phillips  (1959)  would  be  eliminated.  However,  not 
all  terms  in  the  model,  such  as  the  forcing  terms,  are  strictly  linear  or  quadratically  non¬ 
linear.  As  a  result,  some  aliasing  can  occur  in  the  model,  although  results  indicate  that  a 
(3 K  +  1)  x  (3L  +  1)  grid  is  large  enough  to  make  the  non-linear  computational  instability 
negligible. 
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2.3.4  Discretization  in  Time 


The  final  step  in  integrating  the  model  is  the  time  discretization  of  (2.18).  In  the 
present  model,  the  second  order  Adams-Bashforth  scheme  is  used.  We  use  this  particu¬ 
lar  scheme  because  it  damps  the  computational  mode  for  the  oscillation  equation.  This 
property  was  shown  by  Randall  (1989)  and  can  be  proved  as  follows.  Consider  the  linear, 
frictionless  form  of  (2.18),  i.e. 


+  iVklqWklq  =  0. 


(2.32) 


If  we  discretize  (2.32)  using  the  Adams-Bashforth  scheme,  we  can  rewrite  (2.32)  as 

<7"  =  +  »  (|<  -  5<7'))  .  (2.33) 

where  Q  =  ii/ktqAt  and  n  represents  the  time  step.  If  we  assume  (2.33)  has  solutions  of 
the  form  W^+1^  =  AW*]^  then  (2.33)  can  be  written  as 

A2  -  A  (i  +  =  o.  (2.34) 

Equation  (2.35)  has  two  possible  modes  which  can  be  found  through  use  of  the  quadratic 
formula.  These  modes  are 

iSl  +  ^Jl-~U2  +  inj  (2.35) 

and 

in  -  y^i  -  +  .  (2.36) 

Since  A*  — *-  1  as  n  — ►  0,  this  corresponds  to  the  physical  mode.  The  other  mode,  A2, 
approaches  zero  as  fl  — *  0  and  is  therefore  the  computational  mode.  This  is  an  important 
result  because  it  states  that  the  computational  mode  damps  as  A t  becomes  smaller. 

In  order  to  examine  the  behavior  of  the  physical  mode,  it  is  useful  to  simplify  (2.35). 
We  do  this  by  expanding  the  radical  (while  treating  the  terms  containing  f!  within  the 
radical  as  a  single  variable)  using  a  Taylor  series  centered  about  the  point  zero  and  keeping 
the  first  two  terms  in  the  series.  We  are  able  to  do  this  since  is  typically  much  less  than 
one.  The  resulting  expression  can  be  written  as 


Ai  —  1  +  in  ~. 


(2.37) 
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The  modulus  of  (2.37)  is  then 


From  (2.38)  it  is  apparent  that  the  price  paid  for  having  a  damped  computational  mode 
is  the  weak  instability  of  the  physical  mode.  Because  Q  is  raised  to  the  fourth  power  in 
(2.38),  the  instability  of  the  physical  mode  appears  to  be  negligible  if  At  is  sufficiently 
small.  This  assumption  was  verified  through  numerous  model  integrations. 

Before  integrating  (2.18)  it  proves  useful  to  multiply  through  by  the  integrating  factor 
exp(iukiqt).  In  doing  so  (2.18)  can  be  rewritten  as 

jt  [Wu,*****]  =  eiv»'\Fklq  +  Nklq).  (2.39) 

Equation  (2.39)  has  the  property  that  all  linear  terms  may  be  treated  exactly  while  only 
the  nonlinear  and  forcing  terms  need  to  be  treated  explicitly.  Using  this  “semi-exact” 
form,  we  are  able  to  choose  a  time  step  roughly  two  and  a  half  times  greater  than  its 
purely  explicit  counter  part  (DeMaria  and  Schubert,  1984).  If  we  let  t  =  nAt  and  use 
a  forward  time  step  for  the  initial  time  step  and  the  Adams-Bashforth  scheme  for  all 
subsequent  time  steps,  (2.39)  can  be  discretized  as 

WlH  =  WhJ  +  +  F^)]e^At  (2.40) 

and 

=  H'&’e-**""  +  A!  { |  [N$  + 

jWLV1’  +  flf,"  .  (2.41) 

This  is  the  form  of  (2.18)  used  in  the  model  integrations. 

2.4  Advantages  and  Disadvantages  of  Normal  Mode  Spectral  Models 

Perhaps  the  greatest  advantage  of  normal  mode  spectral  models  is  their  ability  to 
partition  the  contribution  to  the  total  flow  by  the  rotational  (slow)  and  gravitational 
(fast)  modes.  For  the  /-plane  model  discussed  above,  this  means  we  can  see  directly 
that  part  of  the  solution  which  is  due  to  gravity-inertia  modes  and  that  which  is  due  to 
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geostrophic  or  balanced  modes.  No  other  type  of  model  can  partition  these  modes  so 
precisely.  In  Chapters  3-5,  we  will  demonstrate  the  usefulness  of  this  feature.  Another 
advantage  of  norms1  mode  spectral  models  is  the  ease  with  which  nonlinear  normal  mode 
initialization  can  be  implemented.  This  results  from  the  models  direct  calculation  of  the 
normal  modes.  Although  this  feature  is  not  used  in  our  simple  model,  it  would  prove  quite 
useful  for  more  complicated  normal  mode  models.  The  final  advantage  of  normal  mode 
spectral  models  is  the  simplicity  with  which  the  time  discretization  can  be  performed.  As 
mentioned  previously,  the  discretization  in  time  reduces  to  a  set  of  ordinary  differential 
equations  for  each  normal  mode  which  can  be  integrated  “semi-exactly.” 

Unfortunately,  the  disadvantages  associated  with  normal  mode  spectral  models  are 
equal  to,  if  not  greater,  than  the  disadvantages.  Perhaps  the  greatest  disadvantage  is 
the  relatively  slow  convergence  rates  of  the  basis  functions  which  must  be  used.  This  is 
especially  true  for  the  cylindrical  coordinate  model  discussed  in  Appendix  A.  In  addition, 
derivatives  of  the  basis  functions  are  not  guaranteed  to  be  eigenfunctions  of  the  linear 
operators  (as  is  the  case  with  the  Cartesian  coordinate  model),  although  they  are  usually 
known  analytically.  This  makes  it  necessary  to  store  large  arrays  containing  these  deriva¬ 
tives  for  each  mode  and  at  each  point  on  the  transform  grid  (see  Appendix  A).  One  more 
disadvantage  is  that  gravity  waves  are  not  allowed  to  simply  radiate  out  the  boundaries 
of  the  model.  They  must  either  be  perfectly  reflected  or,  as  in  the  Cartesian  coordinate 
model,  exit  one  end  of  the  domain  and  enter  the  opposing  end.  Unless  gravity  waves  are 
treated  properly,  they  will  quickly  fill  the  model  domain  with  noise. 

For  the  Cartesian  coordinate  model  which  is  used  in  the  following  chapters,  most 
problems  are  associated  with  the  use  of  a  doubly  periodic  domain.  Special  care  must 
be  taken  to  avoid  the  false  interaction  of  a  disturbance  with  the  same  disturbance  which 
exists  implicitly  in  the  adjacent  domain.  To  help  limit  this  problem,  the  model  domains 
are  kept  larger  than  one  Rossby  radius.  In  linear  dynamics,  this  ensures  that  the  final 
adjusted  geopotential,  associated  with  an  initial  potential  vorticity  (PV)  disturbance,  will 
e-fold  at  least  once  once  before  interacting  with  the  center  of  itself  in  the  adjacent  domain. 
Model  results  indicate  that  this  is  sufficient  to  limit  any  unfavorable  interactions. 


Chapter  3 


INNER  BAND  FORMATION  BY  WAVE  BREAKING  AND  VORTEX 

MERGER 

The  concepts  of  potential  vorticity  (PV)  waves  and  wave  breaking  have  recently  served 
as  useful  tools  in  stratospheric  dynamics  (McIntyre  and  Palmer,  1983,  1984;  Juckes  and 
McIntyre,  1987;  Polvani  and  Plumb,  1992).  Until  now,  however,  these  conceptual  tools 
have  not  spread  to  other  areas  of  meteorology.  In  this  chapter,  we  show  that  these  same 
concepts  also  prove  useful  in  understanding  the  formation  of  a  class  of  hurricane  convective 
bands  which  we  will  refer  to  as  “inner”  bands.  We  begin  by  defining  the  term  “inner”  bands 
and  presenting  observational  evidence  of  their  existence.  We  next  review  the  concept  of 
potential  vorticity  and  mathematically  demonstrate  the  existence  of  PV  waves  in  tropical 
cyclones.  These  concepts  are  then  used  to  understand  the  formation  of  hurricane  bands 
by  the  wave  breaking  process  through  numerical  simulations.  In  addition,  we  analyze  the 
wave  breaking  process  using  the  diagnostic  tools  of  wave-activity  and  wave-activity  flux 
vectors  as  derived  by  Haynes  (1988).  In  the  final  section,  we  discuss  vortex  merger  and 
the  symmetrization  of  PV  anomalies  produced  by  convective  heating. 

3.1  Inner  Bands  Defined 

In  the  present  study,  we  consider  two  types  of  hurricane  convective  bands  -  “inner” 
and  “outer” .  The  latter  will  be  discussed  in  the  following  chapter.  We  define  inner  bands 
to  be  those  convective  bands  which  are  located  typically  less  than  500  km  from  the  storm 
center  and  owe  their  existence  to  asymmetries  in  the  PV  field.  These  asymmetries  can  be 
generated  through  several  processes.  Three  such  processes  are  the  advection  of  the  earth’s 
background  PV,  tropical  cyclone  motion,  and  convection. 
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The  earth’s  background  PV  field  generally  increases  monotonicaUy  with  latitude. 
This  is  due  largely  to  the  earth’s  rotation.  When  the  cyclonic  circulation  associated  with 
a  tropical  cyclone  is  superimposed  on  the  earth’s  natural  PV  field,  the  right  side  of  the 
tropical  cyclone  advects  relatively  low  PV  air  northward  while  the  left  aide  of  the  tropical 
cyclone  advects  relatively  high  PV  air  southward.  The  end  result  is  a  positive  PV  anomaly 
on  the  southern  side  of  the  tropical  cyclone  and  a  negative  PV  anomaly  on  the  northern 
side  of  the  tropical  cyclone.  This  produces  a  wavenumber  one  asymmetry  in  the  PV  field. 
Numerical  simulations  of  this  effect  have  recently  been  presented  by  Shapiro  (1992). 

In  a  similar  fashion,  tropical  cyclone  motion  can  also  induce  a  wavenumber  one  asym¬ 
metry  in  the  PV  field.  As  a  tropical  cyclone  moves  westward  (in  the  northern  hemisphere), 
the  motion  of  the  cyclone  adds  to  the  wind  speed  on  the  northern  side  of  the  vortex  and 
subtracts  from  the  wind  speed  on  the  southern  side  of  the  vortex.  The  vorticity  would 
therefore  be  expected  to  be  slightly  greater  on  the  northern  side  than  on  the  southern 
side.  Since  PV  is  closely  related  to  vorticity,  we  would  again  expect  a  wavenumber  one 
asymmetry  in  the  PV  field. 

In  addition  to  advection  of  the  earth’s  background  PV,  and  tropical  cyclone  motion, 
asymmetries  in  the  PV  field  of  a  tropical  cyclone  can  also  be  generated  by  convection.  If 
we  accept  that  regions  of  convection  are  also  regions  of  high  PV  (as  will  be  discussed  in 
the  following  section),  then  any  asymmetries  in  the  convection  of  a  tropical  cyclone  will 
yield  asymmetries  in  the  PV  field.  An  important  difference  between  this  mechanism  and 
the  previous  two  is  that  convection  is  not  limited  to  wavenumber  one  asymmetries.  All 
wavenumbers  are  possible. 

Asymmetries  in  the  PV  field  can  lead  to  PV  wave  breaking  which  causes  areas  of 
higher  PV  (more  intense  convection)  to  be  ejected  down  gradient  in  the  form  of  a  inner 
bands.  Inner  band  formation  can  also  be  discussed  in  the  context  of  vortex  merger.  That 
is,  when  a  vortex  moves  sufficiently  close  to  a  region  of  relatively  high  PV,  the  vortex  will 
merge  with  the  region  of  high  PV  and  produce  an  inner  spiral  band.  An  example  of  inner 
hurricane  bands  is  presented  in  Fig.  3.1.  The  convective  bands  of  Cleo  (1958)  shown  in 
Fig.  3.1  are  well  within  the  500  km  radius  required  by  our  definition  of  inner  bands.  The 
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Figure  3.1:  An  example  of  inner  hurricane  convective  bands.  Hurricane  Cleo  on  18  August 
1958  when  it  was  just  east  of  Bermuda  and  moving  NNE  at  7  ms-1.  Maximum  winds 
were  approximately  47  ms-1  and  the  minimum  surface  pressure  was  970  mb.  From  Gray, 
1964. 
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processes  of  wave  breaking,  vortex  merger,  and  inner  band  formation  will  be  presented  in 
sections  3.4  and  3.6.  We  will  discuss  the  relation  between  regions  of  high  PV  and  regions 
of  convection  in  section  3.2. 

3.2  Potential  Vorticity  in  the  Shallow  Water  Framework 


In  tropical  cyclone  studies,  whether  observational  or  numerical,  the  output  is  most 
often  examined  in  terms  of  the  wind  and  mass  fields.  In  the  present  study,  we  wish  to 
examine  a  field  which  combines  the  wind  and  mass  fields  to  form  a  single  conservative 
field,  i.e.  PV.  In  this  section,  we  derive  the  PV  equation  in  the  shallow  water  framework 
and  discuss  some  of  its  conservative  properties. 

The  derivation  of  the  PV  equation  for  the  shallow  water  system  is  most  easily  done 
when  the  governing  equations  are  written  in  advective  form.  We  therefore  rewrite  (2.1)- 
(2.3)  as 


du  du  du  ,  d<f>  _ 

(3.1) 

(3.2) 

(3.3) 

where  d/dt  =  d/dt  +  ud/d x  +  vd/dy.  Using  (3.1)— (3.2),  we  first  derive  an  equation  for  the 
vorticity  of  the  fluid.  This  is  done  taking  c?/c?x(3.2)-<?/<?y(3.1).  After  rearranging  terms, 
the  vorticity  equation  can  be  written  as 


, ,  t\  (&**  dv\ 

(/+<)U  +  aj) 


dG_dF 

dx  dy 


(3.4) 


The  PV  equation  can  now  easily  be  obtained  by  combining  (3.3)  and  (3.4)  to  eliminate  the 
divergence  terms.  We  do  this  by  taking  ( H  +  h)(3.4)  -  «  +  /)( 3.3).  After  rearrangement 
of  terms,  the  PV  equation  can  be  written  as 


OF 

dy 


PS\ 
9  ) 


(3.5) 


where  PV  is  represented  by  the  variable  P  which  is  given  as 


P  = 


C  +  f 

H  +  h' 


(3.6) 
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We  see  from  (3.5)  that  PV  is  materially  conserved  for  adiabatic  and  frictionless  flow,  i.e. 
F  =  G  =  5  =  0. 

Another  useful  feature  of  PV  is  that  small  amplitude  gravity  waves  are  invisible  on 
PV-maps.  To  show  this,  we  first  linearize  (3.1)-(3.2)  about  a  basic  state  of  rest  and  derive 
a  PV  equation  in  much  the  same  manner  as  we  did  above.  That  is,  we  first  derive  a  linear 
vorticity  equation  and  eliminate  the  divergence  between  it  and  the  continuity  equation. 
The  resulting  equation  will  be  a  conservation  equation  for  the  perturbation  PV  which  is 
defined  as 

p>  =  c\'  -  ftf.  (3.7) 

One  can  also  show  that  (3.7)  can  be  obtained  (within  a  constant)  by  linearizing  (3.6) 
about  a  basic  state  of  rest. 

We  can  now  easily  show  that  the  gravity  wave  contribution  to  P1  is  identically  zero. 
The  first  step  is  to  transform  (3.7)  to  Fourier  spectral  space  to  obtain 

Hi  =  c2 {ikv'm  -  ilu'M)  -  f4>'k, ,  (3.8) 

where  k,  and  l  again  represent  wavenumbers  in  the  x  and  y  directions  respectively.  The 
second  and  final  step  is  to  refer  back  to  the  transformation  kernel  (2.13)  and  replace  u'kl, 
v'kl  and  <f>'kl  in  (3.8)  with  the  first,  second,  and  third  components  of  K w,  for  q  =  1,2 
respectively.  In  doing  so,  P>  will  sum  identically  to  zero.  We  therefore  conclude  that 
small  amplitude  gravity  waves  are  invisible  on  PV  maps. 

The  PV  discussed  above  is  the  shallow  water  version  of  the  true  Rossby-Ertel  PV 
which  is  defined  as  (1  /p)C,  •  V0  where  p  is  the  density,  <fa  is  the  absolute  vorticity  vector 
and  0  is  the  potential  temperature.  In  regions  of  intense  convection,  where  both  Ca  and  V0 
tend  to  be  large  (especially  their  vertical  components),  the  PV  should  also  be  expected  to 
be  large.  Since  the  shallow  water  PV  is  mathematically  and  physically  analogous  to  the 
Rossby-Ertel  PV,  we  conclude  that  it  too  must  be  large  where  the  convection  is  intense. 
Because  of  this,  we  expect  the  shallow  water  PV  to  not  only  be  a  materially  conservative 
variable,  but  also  an  indicator  of  past  convection.  It  must  be  mentioned  however,  that 
no  computations  of  PV  in  hurricanes  based  on  observed  data  have  been  made  (that  the 
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author  is  aware  of).  This  can  be  explained  by  both  the  lack  of  dense  enough  observational 
data  to  compute  PV  accurately  and  the  lack  of  appreciation  for  PV  in  tropical  dynamics. 

3.3  The  Concept  of  Potential  Vorticity  (Rossby)  Waves  in  Tropical  Cyclones 

The  most  elementary  treatment  of  Rossby  waves  (Platzman,  1968)  begins  with  the 
linearized  (about  a  resting  basic  state)  nondivergent,  barotropic,  /?-plane  model  in  which 
the  northward  gradient  of  basic  state  absolute  vorticity  is  the  constant  (3.  Under  con¬ 
servation  of  absolute  vorticity,  a  sinusoidal  disturbance  causes  Quid  particles  which  are 
displaced  northward  to  acquire  a  clockwise  spin.  A  fluid  particle  which  has  not  been  dis¬ 
placed  and  which  lies  between  northward  displaced  particles  to  its  west  and  southward 
displaced  particles  to  its  east,  will  be  induced  to  move  southward  by  the  clockwise  turning 
particles  to  the  west  and  the  counterclockwise  turning  particles  to  the  east.  Thus,  the 
whole  pattern  will  propagate  to  the  west,  creating  what  is  commonly  referred  to  as  a 
Rossby  wave.  The  mechanism  by  which  Rossby  waves  are  generated  is  referred  to  as  the 
0-effect. 

A  more  general  treatment  of  Rossby  waves  (Hoskins  et  a i,  1985,  section  6)  begins 
with  consideration  of  the  Rossby-Ertel  potential  vorticity  field  displayed  on  an  isentropic 
surface.  A  reasonable  basic  state  is  one  in  which  the  isolines  of  PV  are  oriented  in  the 
east-west  direction  with  higher  values  to  the  north,  which  yields  a  generalized  /7-effect.  If 
the  PV  contours  (which  are  also  material  contours)  are  perturbed  in  a  sinusoidal  fashion, 
a  row  of  alternating  positive  and  negative  PV  anomalies  is  produced.  By  the  invertibility 
principle  (which  can  take  on  a  variety  of  forms  depending  on  the  balance  approximation 
being  used)  these  PV  anomalies  induce  a  flow.  Just  as  in  the  elementary  treatment 
discussed  in  the  previous  paragraph,  the  induced  flow  makes  the  PV  anomalies  propagate 
westward  relative  to  the  basic  flow. 

In  all  but  the  very  upper  tropospheric  part  of  an  axisymmetric  tropical  cyclone  the 
isolines  of  PV  are  circles  with  the  highest  values  of  PV  found  in  the  center  of  the  cyclone 
(Schubert  and  Alworth,  1987).  According  to  the  general  argument  given  in  the  previous 
paragraph,  this  axisymmetric  PV  field  provides  a  basic  state  with  a  monotonic  inward 
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increase  of  PV  on  which  Rossby  waves  can  propagate.  If  the  PV  pattern  is  only  slightly 
disturbed  from  circular  isolines,  there  is  a  restoring  effect  and  Rossby  wave  propagation- 
the  PV  contours  simply  undulate  for  long  periods  of  time.  In  certain  circumstances,  the 
positive  PV  anomalies  quickly  begin  to  elongate  and  form  spiral  tongues  of  high  PV  fluid 
which  become  longer,  thinner  and  more  circular  in  time. 

To  better  understand  PV  waves,  let  us  consider  a  nondivergent,  barotropic  fluid 
system  in  cylindrical  coordinates  which  is  linearized  about  a  basic  state  tangential  flow 
v(r).  In  this  system,  the  vorticity  plays  the  role  of  PV  as  the  conservative  dynamic 
variable.  Following  Lamb  (1945),  the  basic  state  vorticity  is  assumed  to  have  the  top-hat 
form  given  below  in  (3.11).  The  corresponding  basic  state  streamfunction  has  the  form 


^(r)  = 


(  ?£(r2  -  r?)  if  0  <  r  <  ri 
(  ln(r/ri)  if  r%  <  r  <  oo 


(3.9) 


where  £  and  iq  are  constants.  One  differentiation  of  (3.9)  yields  the  basic  state  (Rankine) 
tangential  flow 

di  0  <  r  <  ri 

v(r)  =  ^  ,  ,  (3.10) 

dT  l  hiA/r  if  n  <  r  <  oo 

while  a  second  differentiation  yields  the  piecewise  constant  vorticity  pattern 

d_(  dd}\  [£  if0<r<rx 
rdr 


-  A-  (rAL}  -  lf0<r<r 

rdr  rdr  {  dr)  \Q  ifr,  <r  < 


(3.11) 

ri  <  r  <  oo 

Now  suppose  the  interface  ri  is  perturbed  by  an  amount  r)(<p,  t)  in  the  sinusoidal 
fashion  17(^7,  t)  =  where  17  is  a  complex  constant,  m  the  tangential  wavenumber 

and  v  the  frequency.  The  interface  perturbation  results  in  a  circular  chain  of  vortic¬ 
ity  anomalies,  as  shown  in  Fig.  3.2  for  m  =  4.  Away  from  the  interface  the  vorticity 
is  unperturbed  from  its  basic  state  value  so  that  V2tp'  =  0  for  r  /  rj,  or,  assuming 
^/(r,¥7,<)  =  «(r)e<(m'°-*'*>, 


r—~  -,n2^f=0  for  r^ri+r/.  (3.12) 

dr  \  dr 


As  solutions  of  (3.12)  which  are  1 


'i'(r)  =  A 


{ 


aed  at  r  =  0  and  as  r  — »  oc  we  have 
(r/rj)m  if  0  <  r  <  n  + 17, 

(ri/r)m  if  r\  +  r)  <  r  <  oc, 


(3.13) 
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Figure  3.2:  A  wavenumber  four  perturbation  (dashed)  of  a  circular  (smooth)  PV  surface. 


where  A  is  a  constant  which  must  now  be  related  to  fj.  Note  that  \k(r)  is  a  continuous 
function  in  the  small  amplitude  approximation,  i.e.  as  tj  — *  0.  Continuity  of  v  at  ri  +  rj 
yields 

fcr2  +  mA  c»(mv-*'0  =  i£r2  _  mA 
which,  when  evaluated  at  ri  +  r\  and  linearized,  yields 

mA  =  -  fariii.  (3.14) 


Since  vorticity  is  conserved,  the  vorticity  discontinuity  must  move  with  the  fluid.  For  this 
reason,  the  normal  component  of  velocity  for  a  particle  on  the  boundary  must  be  equal  to 
the  velocity  of  the  boundary  itself.  Using  this,  our  equation  for  the  particle  displacement 
can  be  written 


dri  drt  .  im 


(3.15) 


where  Q  =  vfr.  When  (3.13)-(3.14)  are  substituted  into  (3.15)  we  obtain  the  dispersion 
relation 


(316) 

From  (3.16)  we  see  that  like  Rossby  waves,  these  PV  waves  propagate  against  the  mean 
wind  with  phase  speeds  which  increase  (decrease)  with  wave  length  (number).  Thus, 
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wavenumber  m  =  1  is  stationary  and  wavenumbers  m  =  2,3,4  move  respectively  at  5, 

|  the  speed  of  the  basic  state  tangential  flow.  Since  the  above  analysis  was  first  given  by 
William  Thomson  (Lord  Kelvin)  in  1880  and  was  summarized  by  Lamb  (1932,  pages  230- 
231),  attributing  the  origin  of  this  concept  to  Rossby  (1939)  is  not  historically  accurate. 
However,  calling  the  waves  obeying  the  dispersion  relation  (3.16)  Kelvin  waves  would  be 
confusing  for  obvious  reasons.  Here,  we  shall  follow  conventional  terminology  and  refer  to 
these  waves  as  Rossby  waves,  or  more  generally,  PV  waves. 

3.4  Inner  Band  Formation  by  Wave  Breaking 

The  evolution  of  asymmetric  potential  vorticity  patterns  has  been  studied  extensively 
using  contour  dynamical  methods  for  the  barotropic,  nondivergent  governing  equations 
(e.g.  Dritschel,  1986;  Melander  et  a/.,  1987b).  In  contour  dynamics,  the  evolution  of  a 
two  dimensional  area  of  uniform  PV  (or  vorticity  for  the  case  of  nondivergent,  barotropic 
flow)  is  traced  by  predicting  the  position  of  the  interface  separating  distinct  areas  of  PV. 
This  can  be  done  if  the  flow  is  adiabatic  and  frictionless  since  under  these  conditions  PV 
is  materially  conserved  and  thus  the  PV  boundary  represents  a  material  surface.  In  the 
above  studies,  the  asymmetric  PV  patterns  undergo  a  filamentation  process  which  ejects 
filaments  of  high  PV  air  into  regions  of  relatively  low  PV  air.  This  filamentation  process 
can  be  thought  of  as  an  irreversible  PV  wave  breaking  process  in  which  the  area  where 
the  filaments  exist  is  the  “surf  zone”.  This  terminology  was  first  used  by  McIntyre  and 
Palmer  (1984)  to  describe  the  synoptic  scale  behavior  of  the  winter  middle  stratosphere. 

In  this  section,  we  attempt  to  show  that  the  underlying  dynamics  of  Rossby  wave 
breaking  may  be  used  to  help  understand  the  the  evolution  of  inner  bands  in  tropical 
cyclones.  We  do  this  by  first  discussing  the  concept  of  Rossby  wave  breaking  and  surf 
zones,  and  then  demonstrating  the  usefulness  of  this  concept  through  the  integration  of 
the  shallow  water  model.  The  shallow  water  model  is  useful  because  it  allows  for  max¬ 
imum  resolution  in  the  horizontal  and  maximum  simplification  in  the  vertical  while  still 
capturing  the  fundamental  features  of  vortex  dynamics.  These  features  include  the  nec¬ 
essary  mechanisms  for  PV  waves,  gravity-inertia  waves  and  two-dimensional  turbulence. 


Lastly,  we  demonstrate  that  filamentation  and  wave  breaking  are  slow  manifold  processes, 
thus  suggesting  that  transient  gravity  waves  play  a  very  minor  role  in  the  formation  of 
spiral  bands. 

3.4.1  The  Concept  of  PV  Wave  Breaking  and  Surf  Zones 

As  mentioned  previously,  the  concept  of  wave  breaking  was  first  introduced  observa- 
tionally  by  McIntyre  and  Palmer  (1983,84)  in  their  analysis  of  the  extratropical  middle 
stratosphere.  Numerical  studies  of  this  same  phenomenon  were  later  presented  by  Juckes 
and  McIntyre  (1987).  In  these  studies,  isentropic  surfaces  in  the  stratosphere  were  di¬ 
vided  qualitatively  into  two  distinct  zonally  asymmetric  areas.  These  areas  were  coined 
by  Mclnyre  and  Palmer  as  the  main  vortex,  described  as  an  area  of  high  PV  with  sharp 
PV  gradients  at  its  edge  and  the  surf  zone ,  described  as  an  area  where  systematic  large 
scale  PV  gradients  are  comparatively  weak.  The  name  surf  zone  was  given  to  represent 
the  region  where  PV  wave  breaking  occurs. 

Having  demonstrated  the  existence  of  PV  waves  in  tropical  cyclones  in  the  previous 
section,  we  can  now  define  what  is  meant  by  wave  breaking.  The  concept  of  PV  wave 
breaking  bears  a  conceptual  resemblance  to  waves  which  break  on  the  shores  of  large  bodies 
of  water,  thus  prompting  the  use  of  the  term  surf  zone.  If  a  circular  chain  of  particles  all 
having  the  same  value  of  PV  were  perturbed  in  a  sinusoidal  fashion,  they  would  simply 
undulate  about  their  original  position  and  the  wave  would  propagate  similar  to  a  gravity 
wave  in  a  body  of  water.  Under  certain  conditions,  however,  the  chain  of  particles  does 
not  simply  undulate,  but  rather,  some  particles  sire  ejected  irreversibly  away  from  their 
original  positions  similar  to  waves  breaking  on  the  shore  of  a  beach.  McIntyre  and  Palmer 
(1985)  envisaged  wave  breaking  as  the  rapid,  irreversible  deformation  of  material  contours. 
For  the  adiabatic  and  frictionless  shallow  water  model,  the  material  contours  in  question 
are  simply  isopleths  of  PV.  Wave  breaking  in  this  context  is  important  because  it  often 
determines  whether  PV  waves  will  make  quasi-permanent  changes  to  the  mean  flow  in 
which  they  are  embedded  (McIntyre  and  Palmer,  1985).  Before  continuing,  however,  we 
must  first  define  what  is  implied  by  the  words  “rapid”  and  “irreversible”  in  the  definition 
of  wave  breaking. 
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The  word  “rapid,”  in  the  context  of  wave  breaking,  as  defined  by  McIntyre  and 
Palmer,  implies  the  time  scale  of  the  material  disturbance  should  not  be  much  longer  than 
a  typical  intrinsic  wave  period  in  the  area  of  interest  but  should  also  be  much  shorter  than 
dissipative  time  scales.  We  will  show  later  in  this  section  that  our  interpretation  of  PV 
wave  breaking  in  tropical  cyclones  meets  both  of  these  criteria.  The  word  “irreversible,”  as 
used  by  McIntyre  and  Palmer  is  defined  to  be  conceptually  similar  to  the  definition  used 
in  statistical  mechanics.  That  is,  an  irreversible  process  is  one  in  which  small  changes 
in  the  material  contours  cannot  be  made  to  retrace  their  path.  Unlike  in  statistical  me¬ 
chanics,  however,  we  concern  ourselves  only  with  the  bulk  properties  of  the  fluid,  not  the 
molecular  scale  properties.  As  McIntyre  and  Palmer  point  out,  the  irreversibility  of  the 
material  contours  in  the  fluid  dynamical  sense  would  be  present  regardless  of  whether 
or  not  microscopic  dissipative  terms  such  as  diffusion  were  present.  As  an  example  of 
fluid-dynamical  irreversibility,  consider  Fig.  3.3  which  was  taken  from  Welander  (1955). 
Figure  3.3  shows  the  time  evolution  of  a  passive  material  tracer  in  a  laboratory  tank  con¬ 
taining  a  rotating  fluid.  As  time  moves  forward  from  a-e  the  material  contours  undergo 
an  irreversible  transformation  in  the  sense  that  a  reverse  transformation  from  e-a  is  never 
observed,  even  though  diffusion  plays  only  a  very  minor  role. 

It  should  also  be  mentioned,  as  McIntyre  and  Palmer  (1984)  point  out,  that  wave 
breaking  does  not  imply  instabilities  exist  in  the  flow.  Although  instabilities  can  be  im¬ 
portant  in  certain  cases,  the  irreversible  material  deformations  discussed  above  can  take 
place  in  the  absence  of  any  recognizable  instabilities.  As  we  shall  see  in  chapter  five,  it  is 
the  wave  breaking  process  which  causes  instabilities  and  not  vice  versa. 

Finally,  it  should  be  remembered  that  the  concept  of  wave  breaking  is  not  an  exact 
diagnostic  tool.  There  are  no  mathematical  derivations  or  specific  patterns  which  uniquely 
define  wave  breaking.  The  concept  is  nonetheless  useful  tool  in  understanding  the  complex 
nature  of  wave-mean  flow  interactions.  As  McIntyre  and  Palmer  best  stated, 

“In  reality  we  are  dealing  with  complicated  kinds  of  fluid  motion  to  which 
no  analytical  theory  is  likely  to  be  self-consistently  applicable  over  significant 
spans  of  time,  particularly  linear  theory.  However,  the  lack  of  a  complete 
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Figure  3.3:  An  example  of  fluid-dynamical  irreversibility  produced  by  injecting  dye  into  a 
slowly  rotating  laboratory  tank  of  water.  (From  Welander,  1955.) 
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theory  should  not  deter  one  from  using  a  concept,  especially  one  which  has 
been  found  to  be  useful  as  a  heuristic  organizing  concept  when  trying  to  make 
sense  of  very  complicated  phenomena.”  McIntyre  and  Palmer ,  1985 

To  demonstrate  the  concept  of  wave  breaking,  we  present  a  model  integration  taken 
from  Juckes  and  McIntyre  (1987).  The  model  simulation  was  an  integration  of  the  nondi- 
vergent  barotropic  vorticity  equation  in  spherical  coordinates  with  an  initial  state  consist¬ 
ing  of  a  zonally  averaged  vorticity  field  which  monotonically  increased  with  latitude.  The 
vorticity  field  was  then  forced  with  a  time  dependent  wavenumber  one  disturbance.  The 
purpose  of  the  experiment  was  to  simulate  of  the  effect  of  a  planetary  scale  Rossby  wave 
incident  from  below  on  a  layer  lying  between  two  isentropic  surfaces.  Fig.  3.4  shows  the 
results  from  this  integration.  By  day  four  of  the  integration,  the  wave  breaking  process 
has  begun  to  take  effect  and  the  material  contours  have  begun  a  rapid  and  irreversible 
deformation.  By  day  12  the  main  vortex  has  eroded  considerably  and  the  vorticity  (PV) 
has  spread  over  a  larger  area,  thereby  altering  the  initial  mean  wind  field. 

Although  tropical  cyclone s  are  considerably  different  in  terms  of  scale  and  forcing 
mechanisms  from  the  stratosphere,  the  concept  of  wave  breaking  is  also  of  considerable  use 
in  examining  the  evolution  of  spiral  bands  in  tropical  cyclones.  The  dynamical  similarity 
lies  in  the  fact  that  tropical  cyclones  can  also  be  described  by  a  main  vortex  surrounded 
by  a  surf  zone  in  which  the  spiral  bands  are  found.  In  this  context,  the  wave  breaking 
process  can  be  thought  of  as  modifying  the  mean  wind  profile  by  ejecting  high  PV  air 
down  gradient  into  regions  of  relative  low  PV. 

3.4.2  Wave  Breaking  Experiments 

The  wave  breaking  experiments  presented  in  this  section  are  accomplished  by  specify¬ 
ing  an  asymmetric  initial  vorticity  field  which  is  about  to  undergo  an  irreversible  material 
deformation  and  letting  the  model  evolve  unforced.  We  do  this  to  simulate  the  effects  of 
PV  asymmetries  on  hurricane  band  formation.  In  reality,  hurricane  banding  is  a  contin¬ 
ual  process  in  which  asymmetries  are  constantly  being  generated.  In  contrast,  the  wave 
breaking  experiments  in  this  section  consider  the  evolution  of  a  single  asymmetry  in  the 
PV  field  after  being  formed  by  convective  heating  for  example. 
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Figure  3.4:  An  example  of  the  wave  breaking  process  in  the  mid-stratosphere  as  simulated 
with  a  nondivergent  barotropic  model.  The  initially  symmetric  vortex  is  forced  with  a 
wavenumber  one  disturbance.  Dark  shading  corresponds  to  high  vorticity  (or  PV).  (From 
Juckes  and  McIntyre,  1987.) 
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The  reason  we  choose  to  specify  the  vorticity  field  rather  than  the  PV  is  the  ease  with 
which  it  can  be  inverted  to  obtain  the  initially  balanced  wind  and  mass  fields.  Because 
the  PV  field  closely  mirrors  the  vorticity  field,  we  are  indirectly  specifying  the  shape  of  the 
initial  PV  field  which,  as  mentioned  previously,  is  materially  conserved  for  the  adiabatic 
and  frictionless  flow  considered  in  this  section. 

The  initial  vorticity  patterns  used  in  these  experiments  are  of  two  types.  The  first 
type  of  vorticity  pattern  consists  of  non-intersecting,  non-concentric  circles  whose  radii  all 
fall  along  the  same  axis.  The  second  type  of  vorticity  pattern  consists  of  concentric  ellipses 
all  sharing  similar  orientations  and  aspect  ratios.  These  particular  patterns  were  chosen  to 
represent  primarily  wavenumber  one  and  wavenumber  two  PV  asymmetries  respectively. 
The  methods  by  which  these  fields  are  determined  are  described  below. 

To  obtain  a  wavenumber  one  vorticity  pattern,  we  generate  a  continuous  set  of  noncon- 
centric,  non-intersecting  circles  whose  centers  all  fall  on  the  y-axis.  This  is  a  wavenumber 
one  pattern  in  the  azimuthal  sense;  that  is,  if  a  compass  were  positioned  at  the  center  of 
the  circle  of  smallest  radius,  then  as  the  compass  rotated  360°  one  relative  maximum  and 
one  relative  minimum  would  be  encountered.  In  the  Fourier  sense,  however,  this  would 
only  be  a  primarily  wavenumber  one  pattern.  We  generate  a  vorticity  field  of  this  type 
by  first  defining  the  function  R(r)  which  specifies  how  the  radius  increases  as  the  centers 
of  the  circles  are  moved  away  from  the  position  (xo,  yo)  in  the  center  of  the  domain.  We 
define  this  function  as 

■R(r)  =  —  r  +  y0  (3.17) 

fl 

where  r  is  the  radius  and  rj  is  the  radius  of  the  circle  whose  center  is  located  at  (xo,yi). 
It  should  be  noted  that  if  the  slope  of  (3.17)  exceeds  one,  the  circles  will  intersect.  With 
R(r)  defined,  the  field  of  circles  can  be  defined  using  the  standard  equation  for  a  circle  in 
Cartesian  coordinates,  i.e. 

(x  -  x0)2  +  (y  -  R(r))2  -  r2.  (3.18) 

By  substituting  (3.17)  into  (3.18)  and  using  the  quadratic  formula,  we  can  obtain  the 
radius  of  the  circle  which  satisfies  our  criteria  for  any  point  (x,y).  The  next  step  is  to 
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Figure  3.5:  The  profile  function,  /«(r'),  aa  a  function  of  r'  for  k  =  2.56  (smooth  curve) 
and  k  =  28.38  (dotted  curve). 

specify  the  vorticity,  as  a  function  of  the  radius.  We  follow  Melander  et  af.,  (1987b) 
by  specifying  C(r)  as  a  distribution  with  a  monotone  profile  function  /*.(r),  r  >  0.  We 


therefore  define  the  vorticity  to  be 


(r  <  Ri) 

fK[{r-Ri)nilQ-Ri)}  (ft<r<flo)  , 

(ilo<r) 


(3.19) 


where  Ri  represents  the  radius  at  which  the  vorticity  begins  to  decrease  and  Ro  is  the 
radius  at  which  the  vorticity  vanishes.  The  profile  function,  /«(r),  is  chosen  from  the 
one-parameter  family  {/*;  k  >  0}  where 


/.t(r')  =  exp(-Kr-1  exp(l/(r  -  1))],  0  <  r'  <  1. 


(3.20) 


This  function  smoothly  connects  the  points  Ri  and  Rq.  The  parameter  k  is  the  “steepness 
parameter”  which  determines  how  rapidly  the  vorticity  decreases  as  the  radius  increases. 
If  we  choose  k  =  j(exp2)(ln2)  ~  2.56  then  the  value  of  the  vorticity  will  be  one  half  its 
maximum  value  when  r'  =  0.5.  Likewise,  if  we  chose  k  =  | (exp  4) (In  2)  ~  28.38  then 
the  vorticity  will  be  one  half  its  maximum  value  when  r'  =  0.75.  These  functions  are 
presented  in  Fig.  3.5  for  the  above  mentioned  values  of  k.  We  see  from  Fig.  3.5  that 
region  of  maximum  vorticity  is  much  larger  for  k  =  28.38  and  it  decreases  to  zero  at  a 


much  faster  rate. 
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Using  fche  specified  initial  vorticity  field,  we  next  determine  the  initial  wind  and  mass 
fields  which  are  consistent  with  this  vorticity  field.  We  first  obtain  the  desired  wind  field  by 
assuming  the  wind  is  initially  nondivergent.  For  the  shallow  water  model  this  is  equivalent 
to  assuming  that  dip/dt  =  0.  This  assumption  allows  us  to  use  (3.3)  to  define  a  stream 
function,  ip,  such  that 

u  =  -^,  V  =  ^,  and  c  =  V2t/>.  (3.21)  -  (3.23) 

ay  ox 

By  transforming  £  to  Fourier  spectral  space  the  Laplacian  operator  (3.23)  is  easily  inverted 
to  obtain  the  spectral  coefficients  for  ip  and  its  derivatives,  i.e.  tx  and  v.  Obtaining  the 
physical  space  fields  for  u  and  v  is  then  accomplished  by  simply  performing  an  inverse 
Fourier  transform. 

Having  found  the  appropriate  wind  field,  we  next  determine  the  appropriate  mass 
field.  In  order  for  a  unique  solution  to  exist  and  for  transient  gravity  waves  to  be  mini¬ 
mized,  we  assume  the  wind  and  mass  fields  are  initially  in  balance.  Possibly  the  simplest 
method  of  doing  this  is  through  use  of  the  nonlinear  balance  equation  which  is  obtained 
by  deriving  an  equation  for  the  local  rate  of  change  of  the  divergence  and  setting  it  equal 
to  zero.  A  better,  but  more  complicated,  method  of  obtaining  a  balanced  mass  field  is 
nonlinear  normal  model  initialization.  For  this  method,  the  time  rate  of  change  of  the 
normal  mode  spectral  coefficients  for  q  =  1,2  are  initially  forced  to  zero.  An  excellent 
comparison  of  these  two  methods  can  be  found  in  DeMaria  and  Schubert  (1984).  In  the 
present  study  we  use  the  nonlinear  balance  equation  method. 

The  nonlinear  balance  equation  can  be  derived  by  first  taking  d/dx{3.1)  +  d/dy(2.2). 
By  neglecting  the  local  rate  of  change  of  the  divergence  and  heating  effects  from  the 
resulting  equation,  we  obtain  an  expression  relating  the  wind  and  mass  fields  which  reduces 
the  amount  of  transient  gravity  wave-activity  initially,  i.e. 

=  (3.24) 

Using  (3.24),  the  <p  field  can  easily  be  found  by  inverting  the  Laplacian  in  Fourier  spectral 
space  and  returning  to  physical  space. 
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A  wavenumber  two  vorticity  pattern  is  generated  by  defining  a  continuous  field  of 
concentric  ellipses  with  identical  aspect  ratios  and  orientations.  This  initial  vorticity 
pattern  was  studied  extensively  by  Melander  et  a 1.  (1987b)  with  a  nondivergent  barotropic 
model. 

We  can  easily  generate  an  elliptical  vorticity  field  by  writing  the  equation  for  an  ellipse 
in  slightly  less  familiar  form,  i.e. 

x2  +  ary2  =  r2,  (3.25) 

where  ar  is  the  aspect  ratio  which  is  defined  as  length  of  semi-major  axis  divided  by  the 
length  of  the  semi-minor  axis  and  r  is  now  the  length  of  the  semi-major  axis.  Thus,  if  we 
specify  an  aspect  ratio,  we  can  determine  the  corresponding  length  of  the  semi-major  axis 
for  each  point  (x,y).  Prom  this,  we  can  use  (3.19)— (3.20)  as  we  did  for  the  wavenumber 
one  case  to  determine  the  vorticity  field.  Once  this  is  done  the  wind  and  mass  fields  can 
be  determined  in  the  same  manner  as  discussed  above.  Figure  3.6  shows  the  wavenumber 
one  and  wavenumber  two  type  patterns  generated  by  the  above  formulation  for  k  —  2.56 
(left)  and  k  —  28.38  (right).  The  minimum  contour  is  0.2  with  a  contour  interval  of  0.2. 
In  all  cases,  Rq  —  600  km  and  Ri  =  0  on  a  1600x1600  km  grid.  For  the  wavenumber  one 
patterns,  the  variables  yi  and  ri  were  chosen  to  be  1200  km  and  600  km  respectively. 

3.4.3  Numerical  Results 

The  first  set  of  numerical  integrations  were  performed  using  the  primarily  wavenumber 
one  pattern  described  above  with  no  friction  or  heating.  With  the  absence  of  friction 
and  heating,  PV  is  materially  conserved  and  the  entire  process  becomes  an  advective 
rearrangement  of  PV.  The  purpose  of  this  experiment  is  to  simulate  the  effect  of  an 
asymmetric  tropical  cyclone  in  which  the  convection  is  more  intense  in  the  southern  half 
of  the  domain.  It  needs  to  be  emphasized  that  we  are  in  no  way  assuming  that  heating  and 
friction  are  not  important  aspects  of  a  tropical  cyclone;  they  are  indeed  critical.  We  are 
merely  attempting  to  isolate  a  particular  dynamical  aspect  of  tropical  cyclone  evolution. 

The  model  domain  is  1600  km  by  1600  km  wide  with  2562  gridpoints  on  the  transform 
grid.  In  order  to  evaluate  the  quadratically  nonlinear  terms  exactly,  84  waves  are  kept. 
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Figure  3.6:  Examples  of  wavenumber  one  and  wavenumber  two  patterns  for  different  values 
of  k.  See  text. 
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The  mean  depth,  H,  of  the  fluid  is  300m  giving  a  gravity  wave  speed  (c  =  \fgH)  of 
54.2  ms-1.  The  coriolis  parameter  /  is  chosen  to  be  5.0xl0-5  s_1  which  corresponds  to 
roughly  20°  latitude.  By  taking  c//  we  obtain  a  Rossby  radius  of  deformation  of  1084  km. 
By  choosing  these  values  of  H  and  /  we  insure  that  centers  of  adjacent  vortices  (which 
implicitly  surround  the  vortex  in  question  because  of  the  periodic  boundary  assumption) 
are  at  least  one  Rossby  radius  away  from  each  other.  To  avoid  resolution  problems,  the 
horizontal  scales  of  the  initial  disturbances  are  larger  than  what  would  be  expected  for  a 
typical  tropical  cyclone.  The  fundamental  dynamics,  however,  are  completely  general. 

Figure  3.7  shows  the  initial  wind  and  mass  fields  for  the  first  numerical  integra¬ 
tion  which  has  a  primarily  wavenumber  one  horizontal  structure.  A  maximum  value  of 
5.0xl0~4  s-1  or  10/  was  chosen  for  (  and  Rq  =  300  km  while  Ri  =  0.  The  steepness 
parameter,  k,  was  set  at  2.56  for  this  run.  Note  that  only  the  inner  800x  800  km  part  of 
the  domain  is  shown.  Figures  3.8-3.11  show  the  time  evolution  of  the  fields  at  six  hour 
intervals.  By  six  hours,  wave  breaking  has  clearly  begun.  Wave  breaking  can  be  most 
easily  understood  by  noting  that  the  integration  of  the  vorticity  field  to  get  the  wind 
field  is  a  smoothing  operation.  The  winds  are  therefore  much  more  axisymmetric  than 
the  corresponding  PV  field.  As  a  result,  cross  contour  flow  distorts  the  PV  contours  and 
wave  breaking  takes  place.  As  mass  is  removed  from  the  high  PV  region,  the  PV  gradient 
begins  to  steepen.  This  is  especially  noticeable  in  the  lower  right  quadrant  of  the  domain 
at  6  hrs.  By  twelve  hours  the  region  of  steep  PV  gradient  nearly  encircles  the  main  vor¬ 
tex.  This  process  continues  with  time  until  24  horns,  at  which  point,  the  gradients  have 
become  too  steep  for  the  model  to  resolve  and  the  Gibb’s  phenomenon  begins  to  occur.  It 
is  interesting  to  note  the  remarkable  axisymmetry  of  the  height  field  at  all  times.  Because 
of  the  resolution  of  most  observational  data,  if  a  person  were  to  examine  a  similar  tropical 
cyclone  in  nature,  they  would  undoubtably  miss  the  asymmetries  of  the  PV  field.  This 
suggests  the  need  for  observations  at  high  enough  resolution  that  actual  PV  values  could 
be  accurately  measured. 

Next  we  consider  the  effect  of  wave  breaking  on  the  mean  mean  wind  field.  Shown  in 
Fig.  3.12  are  the  tangentially  averaged  fields  of  vorticity,  normalized  PV,  and  wind  speed 


Figure  3.7:  Dynamic  fields  at  initial  time  for  experiment  one.  Upper  left  panel  corresponds 
to  the  height  field  with  contour  intervals  of  20  m.  Upper  right  panel  corresponds  to  winds. 
The  maximum  wind  vector,  which  is  defined  as  the  distance  between  two  consecutive  tick 
marks,  represents  35  ms-1.  Lower  left  panel  is  the  normalized  PV  (i.e.  P  x  H).  Lower 
right  panel  is  the  absolute  vorticity  field.  The  normalized  PV  and  absolute  vorticity  fields 
are  in  units  of  10~5  s'1  and  have  contour  intervals  of  l.OxlO-4  s'1. 
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Figure  3.8:  Same  as  Fig.  3.7  but  for  6  hours. 
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Figure  3.9:  Same  as  Fig.  3.7  but  for  12  hours. 


Figure  3.10:  Same  as  Fig.  3.7  but  for  18  hours. 
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Figure  3.11:  Same  as  Fig.  3.7  but  for  24  hours. 
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Figure  3.12:  Tangentially  averaged  values  of  vorticity,  normalized  PV,  and  wind  speed  at 
18  hours  for  experiment  one.  Dotted  curves  correspond  to  initial  time  and  smooth  curves 
correspond  to  18  hours.  The  normalized  PV  and  vorticity  are  in  units  of  10-6  s-1  and 
the  wind  speed  is  in  ms-1. 
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at  the  initial  time  (dotted  curves)  and  at  18  hours  (smooth  curves).  Although  the  PV 
is  redistributed,  the  effect  of  wave  breaking  on  the  tangential  mean  is  hardly  noticeable. 
This  is  also  true  for  the  tangentially  averaged  wind  and  vorticity  fields.  In  the  following 
wave  number  two  experiment,  we  will  see  that  this  is  not  always  the  case. 

In  the  second  experiment,  the  model  is  initialized  with  a  primarily  wavenumber  two 
vorticity  pattern.  As  with  the  previous  experiments,  the  model  run  is  to  isolate  the 
fundamental  dynamics  of  a  tropical  cyclone  in  which  the  convective  heating  has  a  primarily 
wavenumber  two  horizontal  structure.  The  initial  vorticity  for  experiment  two  has  a 
maximum  vorticity.  £,  of  6.0 x  10-4  s-1  or  12/.  Other  parameter  settings  include:  ar  — 
2.0,  k  =  2.56,  Ro  ss  400km,  and  R,  =  0.  This  vorticity  field  corresponds  to  a  vortex 
considerably  larger  than  that  expected  for  a  typical  hurricane;  however,  the  results  are 
completely  general  for  the  simple  model  we  are  considering  here.  The  reason  for  the  large 
vortex  is  to  maximize  the  horizontal  resolution  while  insuring  the  domain  is  at  least  one 
Rossby  radius  in  length.  The  initial  fields  for  experiment  two  are  shown  in  Fig.  3.13.  From 
Fig.  3.13  we  again  see  that  the  height  field  is  considerably  more  symmetric  than  the  PV 
field.  This  again  sets  the  stage  for  wave  breaking.  Figures  3.14-3.17  show  the  evolution 
of  the  fields  for  a  24  hour  period.  By  six  hours  wave  breaking  has  clearly  begun.  In 
contrast  to  experiment  one,  the  wave  breaking  process  causes  two  tongues  of  PV  to  be 
ejected,  one  at  each  end  of  the  semi-major  axis.  As  time  proceeds,  these  tongues  elongate 
forming  what  we  refer  to  as  inner  bands.  By  18  hours  the  steepening  of  the  PV  gradient 
by  the  wave  breaking  process  has  caused  Gibb’s  phenomena  to  again  occur.  Since  no 
damping  terms  were  included  in  the  integration,  the  effect  of  this  is  to  generate  spurious 
noise  in  the  calculations  beyond  this  point.  By  24  hours,  the  bands  completely  surround 
the  the  main  vortex. 

Another  prominent  feature  of  experiment  two  is  the  symmetrization  of  the  initially  el¬ 
liptical  PV  pattern.  These  results  are  consistent  with  the  nondivergent  results  of  Melander 
et  a 1.  (1987b).  Heuristically,  the  symmetrization  process  can  be  understood  by  examining 
the  radial  wind  field.  Figure  3.18  shows  the  initial  radial  wind  field  superimposed  over  the 
initial  PV  field.  As  the  tangential  wind  acts  to  rotate  the  PV  pattern  cyclonically,  it  moves 


Figure  3.13:  Initial  fields  for  experiment  two.  Panels  and  contouring  are  the  same  as  for 
experiment  one. 
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Figure  3.14:  Same  as  Fig.  3.13  but  at  6  hours. 


52 


Figure  3.15:  Same  as  Fig.  3.13  but  at  12  hours. 
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Figure  3.16:  Same  as  Fig.  3.13  but  at  18  hours. 
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Figure  3.17:  Same  as  Fig.  3.13  but  at  24  hours. 
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Figure  3.18:  The  radial  wind  and  PV  fields  at  the  initial  time  for  experiment  two. 

the  regions  relatively  weak  PV  gradient  (either  end  of  the  semi-major  axis)  into  areas  of 
inward  radial  flow.  Likewise,  regions  of  relatively  sharp  PV  gradients  (either  end  of  the 
semi-minor  axis)  are  being  moved  into  areas  of  outward  radial  flow.  The  end  effect  is  to 
tighten  the  gradient  where  it  is  weak  and  loosen  it  where  its  strong,  thus  symmetrizing 
the  PV  pattern. 

We  now  consider  the  effects  of  the  wave  breaking  and  symmetrization  processes  on 
the  maintenance  of  the  mean  winds.  Figure  3.19  shows  the  tangentially  averaged  vorticity, 
normalized  PV  and  wind  speed  at  the  initial  time  (dotted)  and  at  18  hours  (smooth)  for 
experiment  two.  The  symmetrization  process  has  caused  the  tangentially  averaged  PV 
to  slightly  increase  within  180  km  of  the  vortex  center.  In  addition,  the  symmetrization 
process  has  caused  the  tangentially  averaged  PV  to  decrease  slightly  between  180  km  and 
300  km  radius.  Beyond  300  km,  although  it  is  not  clearly  evident  from  the  figure,  the 
tangentially  averaged  PV  has  again  increased  due  to  ejection  of  higher  PV  air  to  large 
radii  by  the  wave  breaking  process.  Although  there  is  no  direct  mathematical  correlation 
between  PV  and  wind  speed,  the  trend  in  the  tangentially  averaged  PV  field  closely  mirrors 
the  trend  in  tangentially  averaged  vorticity  (except  at  small  radii)  with  which  there  is  a 
direct  mathematical  correlation,  i.e.,  (  =  d(rv)/rdr.  As  the  PV  and  vorticity  are  shifted 
outwards,  there  is  both  an  increase  in  the  maximum  wind  speed  as  well  as  an  outward 


Figure  3.19:  Tangentially  averaged  vorticity  (top),  normalized  PV  (middle)  and  wind 
speed  at  initial  time  (dotted)  and  18  hours  (smooth).  The  normalized  PV  and  vorticity 
are  in  units  of  10~6  s-1  and  the  wind  speed  is  given  in  ms-1. 
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shift  in  the  radius  of  maximum  winds.  Correspondingly,  the  winds  decrease  at  large  radii. 
Dynamically,  we  should  expect  higher  tangentially  averaged  wind  speeds  when  the  area 
under  the  tangentially  averaged  vorticity  curve  is  higher.  This  can  easily  be  understood 
by  integrating  the  above  vorticity  relation  for  Rankine  vorticity  distributions  which  are 
equally  wide  but  differing  in  maximum  vorticity.  The  vorticity  and  wind  speed  curves 
shown  in  Fig.  3.18  also  demonstrate  this  relationship.  Moving  radially  outward  from 
the  origin,  the  area  under  the  initial  vorticity  field  is  slightly  larger  than  the  area  under 
the  vorticity  curve  at  18  hours;  correspondingly,  the  winds  are  slightly  larger  initially. 
However,  at  ~110  km  the  area  under  the  vorticity  curve  at  18  hours  begins  to  exceed 
the  area  under  the  initial  vorticity  curve  and  the  winds  at  18  hours  are  correspondingly 
larger.  At  ~240  km  the  pattern  is  again  reversed  as  the  area  becomes  larger  under  the 
initial  vorticity  curve. 

In  the  above  discussions  of  both  wave  breaking  and  axisymmetrization,  we  treated 
PV  as  if  it  were  a  passive  material  tracer.  We  can  easily  demonstrate  that  this  is  true 
by  simulating  the  evolution  of  a  passive  material  tracer  which  has  the  same  initial  profile 
as  the  PV  field.  This  is  accomplished  simply  by  integrating  dTr/dt  —  0,  where  Tr  is  any 
material  tracer,  using  the  model  winds  to  advect  the  tracer.  Results  from  this  test  (not 
shown)  show  that  the  tracer  evolves  in  an  identical  fashion  to  the  PV.  This  should  not 
be  surprising,  since  PV  was  shown  to  be  materially  conserved  in  section  3.2.  The  fact 
that  PV  behaves  simply  as  a  passive  material  tracer  in  the  adiabatic,  frictionless  case, 
reemphasizes  its  usefulness  as  a  dynamic  variable. 

Passive  tracer  experiments  are  also  useful  for  viewing  the  Lagrangian  evolution  of 
spiral  bands.  That  is,  they  allow  us  to  view  particle  trajectories  which  in  turn  allows  us 
to  determine  where  the  air  in  a  spiral  band  originated.  We  do  this  by  defining  two  passive 
tracer  fields,  one  which  initially  varies  as  sin[(27r/£I)x]  and  one  which  initially  varies  as 
sin((2ir/Irv)y].  When  displayed  graphically,  these  two  tracers  create  sets  of  perpendicular 
lines,  respectively,  which  range  from  -1  to  1  over  our  inner  800x800  km  display  region. 
Clearly,  fields  which  provide  equally  spaced  contours  would  be  preferred,  however,  we  chose 
sine  functions  because  they  satisfy  the  models  periodic  boundary  conditions.  By  displaying 
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these  two  fields  together,  intersecting  isopleths  uniquely  define  a  single  particle  which  can 
easily  be  tracked  in  time  by  noting  the  position  of  the  intersection.  The  evolution  of  the 
passive  tracer  field  is  completely  decoupled  horn  the  dynamics  of  the  system;  that  is,  the 
model  winds  are  used  only  to  provide  the  advection.  Figure  3.20  shows  the  deformation 
of  this  passive  tracer  field  by  the  fiow  field  associated  with  experiment  two.  In  addition  to 
material  conservation,  the  total  mass  bounded  by  any  four  isopleths  is  also  be  conserved. 
We  have  therefore  shaded  five  initial  squares,  two  which  are  initially  located  on  the  lobes 
of  the  ellipse  (i.e.  either  end  of  the  semi-major  axis),  two  which  are  initially  located  off 
the  ellipse,  and  one  in  the  center.  As  time  evolves,  the  passive  tracer  field  is  stretched 
as  it  wraps  around  the  main  vortex.  Although  the  four  uncentered  regions  are  initially 
located  equal  distances  from  the  center  of  the  vortex,  the  two  which  were  initially  located 
on  the  semi-major  axis  are  quickly  stretched  and  wrapped  around  the  main  vortex  while 
the  two  which  were  initially  located  off  the  ellipse  are  stretched  to  a  much  lesser  extent  and 
remain  separated  by  nearly  equal  distances.  It  appears  as  if  the  symmetrization  process 
is  accomplished  by  advecting  the  some  of  the  high  PV  from  the  lobes  of  the  ellipse  into 
the  center  of  the  vortex.  Squares  located  toward  the  center  of  the  vortex,  however,  simply 
rotate  with  very  little  distortion.  This  is  not  surprising  since  the  vorticity  near  the  center 
is  very  close  to  that  of  solid  body  rotation.  It  is  also  interesting  to  note  the  evolution  of 
the  dot  Fig.  3.20  which  represents  a  single  fluid  particle.  This  particle  is  initially  located 
near  the  right  end  of  the  semi-major  axis.  By  day  six,  this  particle  is  located  near  the 
tip  of  the  incipient  band.  This  allows  us  to  say  for  certain  that  the  air  in  both  bands 
originated  near  the  lobes  of  the  ellipse. 

Lastly,  we  consider  the  contributions  to  the  total  solution  by  the  slow  and  fast  modes. 
This  is  done  simply  by  summing  over  either  the  q  =  0  coefficients  (slow  modes)  or  the 
q  =  1,2  coefficients  (fast  modes)  as  mentioned  in  chapter  two.  The  effect  of  partitioning 
the  flow  is  shown  in  Figs.  3.21-3.22  which  give  the  fast  and  slow  mode  contributions  to 
the  vorticity  field  (Fig.  3.21)  as  well  as  the  wind  and  mass  field  (Fig.  3.22).  It  should 
be  mentioned  that  the  partitioning  of  nonlinear  quantities  such  as  PV  is  also  possible, 
however,  the  sum  of  the  slow  and  fast  mode  contributions  will  not  give  the  total  solution. 


Figure  3.20:  The  evolution  of  a  material  tracer  field  for  experiment  two  at  initial  time 
(upper  left),  2  hours  (upper  right),  4  hours  (lower  left)  and  6  hours  (lower  right).  Contours 
are  from  -0.9  to  0.9  by  0.2.  Shaded  areas  are  to  aid  in  viewing  particle  trajectories. 


Figure  3.21:  Slow  and  fast  mode  contributions  to  the  vorticity  field  at  12  hours  for  exper¬ 
iment  one.  Total  solution  is  shown  on  the  top  while  the  slow  and  fast  mode  contributions 
are  shown  on  the  lower  left  and  right  respectively.  Contours  are  the  same  as  in  Fig.  3.7 


Figure  3.22:  Slow  and  fast  mode  contributions  to  the  wind  and  mass  fields  at  12  hours 
for  experiment  one.  Total  solution  is  shown  on  the  top  while  the  slow  and  fast  mode 
contributions  are  shown  on  the  lower  left  and  right  respectively.  Contours  are  the  same 
as  in  Fig.  3.7 
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Figure  3.23:  A  schematic  diagram  of  Leith’s  slow  manifold.  See  text  for  explanation. 

For  this  reason,  we  do  not  partition  nonlinear  variables.  It  clear  from  Fig.  3.21  that  gravity 
waves  play  a  very  minor  role  in  the  banding  process.  This  is  not  to  say  that  hurricane 
bands  are  completely  balanced  phenomena  either.  This  fact  is  evident  from  the  extensive 
analysis  of  Ryan  et  a I.  (1992)  and  Powell  (1990).  In  their  analysis,  Ryan  et  a J.  found 
small  deviations  of  the  D  value  associated  with  the  bands.  The  D  value  is  defined  as 
the  difference  between  the  radar  altitude  and  the  pressure  altitude  with  reference  to  the 
standard  atmosphere.  It  is  therefore  analogous  to  the  deviation  in  geopotential  from  its 
mean  value.  The  band  was  found  to  be  an  area  of  increased  low  level  convergence  (see 
Chapter  1).  Similar  results  were  also  shown  by  Powell  (1990).  Both  of  these  findings 
suggest  the  presence  of  gravity  wave  modes.  In  addition,  Ryan  et  a /  and  Powell  both  show 
the  bands  to  be  associated  with  high  relative  vorticity,  just  as  our  model  predicted.  This 
is  suggestive  of  rotational  modes.  It  therefore  seems  inappropriate  to  classify  hurricane 
bands  as  either  purely  rotational  or  purely  gravity  wave  phenomenon.  A  more  appropriate 
description  is  that  hurricane  bands  are  slow  manifold  phenomena,  as  defined  by  Leith 
(1980),  which  project  onto  both  the  slow  and  fast  modes. 

The  lack  of  deviations  from  axisymmetry  in  the  height  (or  geopotential)  field  in  our 
model  simulations  indicates  that  the  bands  created  in  the  model  are  not  truly  representa¬ 
tive  of  those  in  nature.  One  explanation  for  this  difference  is  that  the  model  evolves  near 
the  adiabatic  and  frictionless  slow  manifold,  whereas  in  nature,  hurricanes  evolve  along 
the  diabatic  and  frictional  slow  manifold.  Figure  3.23  shows  schematically  both  the  adia¬ 
batic  and  frictionless  slow  manifold  as  well  as  the  diabatic  and  frictional  slow  manifolds. 
Drawn  on  the  abscissa  is  the  two  dimensional  representation  of  the  multi-dimensional  ro- 
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tational  manifold  (R)  which  is  the  set  of  all  linearly  independent  eigenvectors  that  can  be 
combined  to  form  any  purely  rotational  field.  Likewise,  the  ordinate  is  a  two-dimensional 
representation  of  the  multi-dimensional  gravitational  manifold  (G)  which  is  the  set  of  all 
linearly  independent  eigenvectors  that  can  be  combined  to  form  any  purely  gravitational 
field.  The  slow  manifold,  whether  adiabatic  and  frictionless  or  diabatic  and  frictional,  is 
defined  (to  a  low  order  approximation)  to  be  the  locus  of  all  points  such  that  the  total 
derivative  of  all  gravitational  modes  is  zero.  Thus,  evolution  along  the  slow  manifold  is 
characterized  by  the  absence  of  transient  gravity-wave  activity.  In  this  sense,  the  gravity- 
wave  modes  are  enslaved  to  the  balanced  modes.  By  using  the  nonlinear  balance  equation 
to  initialize  our  experiments,  we  are  allowing  our  model  to  evolve  near  the  adiabatic  and 
frictionless  slow  manifold.  To  evolve  exactly  on  the  slow  manifold  would  require  the  use  of 
nonlinear  normal  mode  initialization  (Leith,  1980;  DeMaria  and  Schubert,  1984).  Nature, 
however,  prefers  to  evolve  near  the  diabatic  and  frictional  slow  manifold  in  which  gravity- 
wave  activity  is  more  prominent.  The  reason  for  the  increased  gravity-wave  activity  is  the 
addition  of  frictional  convergence  and  diabatically  driven  convection.  This  gravity-wave 
activity  would  still  be  expected  to  be  enslaved  by  the  rotational  modes;  that  is,  transient 
gravity-waves  are  expected  to  be  minimal.  In  short,  spiral  bands  project  onto  both  the 
gravitational  and  rotational  modes,  thus  pure  gravity  wave  theories  are  not  completely 
correct.  It  is  likely  that  evolution  along  the  diabatic  frictional  slow  manifold  provides  the 
observed  deviations  from  axisymmetry  in  the  geopotential  field.  In  the  final  section  of 
this  chapter,  we  will  show  that  the  addition  of  a  diabatic  heat  source  produces  bands  and 
geopotential  fields  which  are  more  representative  of  nature. 

Another  reason  for  partitioning  the  wind  and  mass  fields  is  to  examine  roles  the 
slow  and  fast  modes  play  in  providing  gradient  balance.  As  pointed  out  by  Schubert 
and  DeMaria  (1985)  using  am  axisymmetric  model,  the  gravity  wave  contribution  to  the 
total  solution  yields  ur  «  0,  «  0,  and  0/0  and  is  responsible  for  the  cyclostrophic 

part  of  gradient  balance.  This  is  also  clearly  observed  in  Fig.  3.22.  In  terms  of  the  slow 
manifold  diagram  (Fig.  3.23),  the  gravity-wave  contribution  takes  us  from  a  point  on 
the  rotational  manifold  to  a  point  on  the  adiabatic  and  frictionless  slow  manifold.  The 
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addition  of  diabatic  and  frictional  forcing  would  then  move  us  to  a  point  on  the  diabatic 
and  frictional  slow  manifold. 

3.5  Wave- Activity  and  Wave-Activity  Flux  Diagnostics 

In  the  previous  section,  we  examined  the  evolution  of  inner  bands  in  terms  of  the 
dependent  variables  as  well  as  PV.  We  now  consider  the  use  of  wave-activity  as  a  locally 
conserved  variable. 

If  the  atmosphere  is  partitioned  into  mean  and  eddy  components,  it  can  often  be 
shown  that  the  eddy  quantities  will  satisfy  a  local  conservation  relation  of  the  form 

dA 

^  +  VF  =  5,  (3.26) 

where  both  A  and  F  are  functions  of  eddy  quantities.  The  variable  A  is  the  wave-activity 
and  F  represents  a  flux  of  the  wave-activity.  The  right-hand  term,  5,  is  then  a  source  or 
sink  of  A.  The  first  to  derive  an  expression  similar  to  (3.26)  were  Eliassen  and  Palm  (1960) 
for  steady  conservative  waves  (i.e.  dA/dt  —  S  =  0)  in  the  quasi-geostrophic  framework. 
It  is  for  this  reason  that  the  flux  term  in  (3.26)  is  commonly  referred  to  as  the  Eliassen- 
Palm  (EP)  flux.  To  distinguish  between  the  steady  conservative  waves  considered  by 
Eliassen  and  Palm  and  more  general  applications,  we  will  refer  to  F  simply  as  the  flux 
of  wave-activity.  For  quasi-geostrophic  flows,  the  flux  of  wave-activity  can  be  directly 
related  to  the  zonal  wind.  It  therefore  provides  an  excellent  tool  for  the  study  of  wave- 
mean  flow  interactions.  Andrews  and  McIntyre  (1978a, b)  were  able  to  generalize  these 
concepts  to  the  primitive  equations,  but  their  formulation  required  knowledge  of  particle 
displacements,  thus  making  the  theory  difficult  to  use  in  practice.  McIntyre  and  Shepherd 
(1987)  eliminated  the  need  for  particle  displacements  and  derived  conservation  equations 
above  for  the  nondivergent  barotropic  equations  which  were  strictly  in  terms  of  Eulerian 
quantities.  In  a  recent  work,  Haynes  (1988)  extended  the  finite  amplitude  wave  activity 
conservation  laws  of  McIntyre  and  Shepard  to  the  hydrostatic  primitive  equations.  The 
theory  is  incomplete,  however,  in  the  sense  that  the  flux  of  wave-activity  cannot  be  directly 
related  to  changes  in  the  mean  flow.  It  is  nonetheless  a  useful  quantity  to  consider  because 
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of  its  conservative  properties  which  are  directly  related  to  wave  amplitude.  In  addition, 
this  is  a  first  step  towards  understanding  wave-mean  flow  interactions. 

In  this  section,  we  simplify  Haynes’  results  to  the  unforced  shallow  water  equations. 
These  equations  are  then  used  to  calculate  the  wave-activity  and  the  associated  EP  flux 
vectors.  We  begin  by  an  in-depth  derivation  of  the  wave-activity  relation  using  the  mo¬ 
mentum  Casimir.  The  derivation  is  followed  by  numerical  calculations. 

3.5.1  The  Wave-Activity  Relation  for  the  Shallow  Water  Equations 


We  consider  a  shallow  water  system  on  an  /-plane.  Using  cylindrical  coordinates 
in  the  horizontal  and  denoting  the  radial  wind  by  u  and  the  tangential  wind  by  v,  the 


governing  set  takes  the  form 


d(hu)  d(rhuu)  d(kvu) 


~{f+fjhv  +  gh^  =  0,  (3.27) 


d[h(rv  +  \fr2)]  ,  d[rhu(rv  +  \fr2)\  d[hv(rv  +  ^fr2)]  ,  Ldh  n 

- dT - + - 7dT - + - - +9%  =  °’  (328) 

where  h  now  represents  the  fluid  depth.  We  regard  (3.27)-(3.29)  as  a  closed  system  in 
u,  v  and  h.  A  consequence  of  (3.27)-(3.29)  is  the  conservation  of  the  potential  vorticity 
P  ~  (C  4-  f)/h.  Since  P  is  conserved,  any  function  C(P)  is  also  conserved.  Thus,  the 


Casimir  function  C(P)  satisfies 

d(hC)  djruhC)  d(vhC) 

dt  rdr  "**  rd<p 

Adding  (3.28)  and  (3.29),  we  obtain 

d[h(rv  +  j/r2  +  C)]  d[ruh(rv  +  5/r2  +  C)]  d[vh(rv  +  ^fr2 +  C) +  r^gh2] 
dt  +  rdr  +  rd<p 


(3.30) 


(3.31) 


Following  McIntyre  and  Shepherd  (1987)  and  Haynes  (1988),  we  now  wish  to  choose 
C{P)  such  that  the  disturbance  part  of  h(rv  +  5/r2  +  C)  can  be  written  as  a  divergence, 
plus  a  term  which  is  explicitly  second  order  in  wave  amplitude.  The  disturbance  part  of 
h(rv  +  5/r2  +  C)  is  defined  as 

{h  [rv  4-  5/r2  +  C{P)  1  =  h  ^rv  4-  \fr2  4-  C(P)J  -  ho  [r«o  4-  5/r2  4-  C(Po)]  • 
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After  some  manipulation  this  can  be  written 

{h  [rv  +  \fr 2  +  C(P)]  }e  =  rvehe  +  h  (C(P)  -  C(P0)  -  PeC'(Po)] 

d[C(P0)rve]  a[C(P0)tie] 
rdr  rd<p 

+he  [r«o  +  £/r2  +  C(P0)  -  PqC (Po)]  +  rve  [ho  -  •  (3-32) 

The  right  hand  side  of  the  first  line  of  (3.32)  is  second  order  in  disturbance  amplitude  and 
is  defined  to  be  the  wave-activity.  The  second  line  is  the  divergence  of  a  flux.  We  now 
wish  to  choose  C(P)  in  such  a  way  that  the  last  line  in  (3.32)  vanishes.  The  choice 

fP 'mm* 

C(P)  =  -Jp  m(P)dP,  (3.33) 

where 

m(Po(r))  =  f  ho(f)fdf,  (3.34) 

Jo 

makes  the  last  term  in  the  square  brackets  in  (3.32)  vanish,  as  can  easily  be  verified.  Since 

^  [rvo  +  \fr2  +  C(P0)  -  PoC^Po)]  *  Po  [ho  -  ,  (3-35) 

and  since  the  right  hand  side  of  (3.35)  vanishes  by  the  previous  argument,  we  conclude 
that  rv 0  +  5/r2  +  C(Po)  —  PqC(Po)  is  constant  with  radius.  The  constant  is  zero  because 
the  left  hand  side  vanishes  at  r  =  0.  Thus,  the  entire  third  line  in  (3.32)  vanishes  and  we 
can  write  (3.32)  as 

{A  [r»  +  +  0(J»)] }.  =  .4  +  222^1^1  _  (3.36) 

where 

A  =  rvehe  +  h  [C(P)  -  C(P0)  -  Pem(P0)] 

is  the  wave-activity. 

Noting  that  the  basic  state  is  a  steady  solution  of  the  governing  equations  and  then 
substituting  (3.36)  into  (3.31)  we  obtain 


tr + 7k  ['ir“(r” + + c> + 
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H™ + 5/r2 +c)+ ■ m(po)^] = °- 

Using  the  disturbance  momentum  equations 


%  -  hPv + wvo + = 0, 
at  or 

d(™e)  +  ^  +  d(ghe  +  %uj  +  VQVe  +  \vp  _  0 
dt  d<p 

in  (3.37),  we  obtain 

dA  d[r(uA  +  hQrueve)\  d[vA  +  hor^jv?  -  u?)  +  r\ghl\ 
dt  rdr  rd<p 

where  the  wave-activity  can  also  be  written  as 


A  =  rvehe  +  [m(P)  ~  m(P0)]  dP. 


(3.37) 

(3.38) 

(3.39) 

(3.40) 

(3.41) 


3.5.2  Calculation  of  Wave- Activity  and  Wave- Activity  Flux  Vectors 

Before  presenting  calculations  of  wave-activity  and  its  flux,  it  is  important  to  discuss 
the  usefulness  of  such  quantities.  Unlike  generalized  Lagrangian  mean  theory,  the  wave- 
activity  flux  cannot  be  directly  correlated  to  changes  in  the  mean  flow,  nor  can  we  make 
statements  about  the  positive  definiteness  or  negative  definiteness  of  A.  What  we  can 
say,  however,  is  that  wave-activity  offers  a  calculable  and  conservative  measure  of  wave 
amplitude.  In  the  absence  of  diabatic  or  frictional  effects,  this  implies  that  convergence 
of  the  flux  vectors  cannot  occur  unless  the  wave  amplitude  is  building  up  transiently  in 
the  area  of  convergence.  For  our  present  purposes,  we  should  expect  banded  areas  to  also 
be  areas  of  high  wave-activity.  In  addition,  convergence  of  the  wave-activity  flux  vectors 
should  correspond  to  regions  of  increasing  wave-activity. 

To  demonstrate  the  usefulness  of  wave-activity,  we  have  calculated  both  A  and  its 
flux  vectors  as  defined  in  (3.40)  for  the  wave-number  two  experiment  discussed  previously. 
Again  the  model  is  integrated  using  2562  points  on  the  transform  grid  and  keeping  84 
waves.  The  calculations  for  A  are  performed  using  the  definition  given  in  (3.36).  The 
mean  state  was  computed  as  the  tangential  average  at  the  initial  time.  These  averages 
were  computed  on  a  cylindrical  grid  with  a  maximum  radius  of  600  km.  The  fields  were 
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sampled  every  12  km  in  the  radial  direction  (60  grid  points)  and  every  6°  in  the  tangential 
direction  (60  grid  points).  Since  the  tangentially  averaged  initial  fields  are  axisymmetric 
and  in  gradient  balance,  they  provide  a  steady  solution  of  the  governing  equations.  It 
should  be  noted  that  because  the  calculations  are  done  in  cylindrical  coordinates,  the 
plotting  of  data  requires  the  interpolation  of  the  data  back  to  Cartesian  coordinates. 
Because  of  this,  the  plots  appear  to  be  slightly  less  smooth  than  previous  figures.  Figure 
3.24  shows  both  the  wave-activity  and  the  flux  vectors  for  a  12  hour  integration.  Clearly, 
the  inner  bands  corresn  >nd  to  regions  of  high  wave-activity. 

Initially,  the  wave-activity  is  the  highest  where  the  deviation  of  the  PV  field  from 
axisymmetry  is  the  greatest.  From  (3.40),  we  see  that  A  satisfies  a  local  conservation 
relation.  This  implies,  the  area  integral  of  A  is  constant.  An  interesting  observation  of 
the  wave-activity  is  that  not  only  does  it  satisfy  a  local  conservation  relation,  but  it  also 
appears  to  be  nearly  materially  conserved.  That  is,  it  appears  to  follow  the  same  trend 
as  the  material  tracer  trajectories  of  the  previous  section. 

Since  the  wave-activity  relation  used  in  here  is  not  part  of  a  complete  theory,  it 
cannot  be  said  for  certain  that  there  is  a  relationship  between  the  wave-activity  flux  and 
wave-mean  flow  interactions. 

3.6  Inner  Band  Formation  by  Vortex  Merger 

In  the  previous  sections  we  considered  the  formation  and  evolution  of  inner  bands 
through  wave  breaking  processes.  We  now  demonstrate  how  spiral  bands  can  also  be 
formed  when  a  pre-existing  cyclonic  vortex  merges  with  a  region  of  relatively  high  PV.  Vor¬ 
tex  merger  is  not  a  new  concept  and  has  been  examined  in  numerous  fluid  dynamical  stud¬ 
ies  (e.g.  McWilliams,  1984;  Griffiths  and  Hopfinger,  1987;  Melander  et  aJ.,  1987a, 1988). 
The  application  of  vortex  merger  to  tropical  cyclones,  however,  is  new.  Recent  work  in 
this  area  has  been  accomplished  by  Evans  and  Holland  (1991),  Holland  and  Evans  (1992), 
Holland  and  Lander  (1992),  Lander  and  Holland  (1992),  Ritchie  and  Holland  (1992),  and 
Holland  and  Dietachmayer  (1992).  In  the  present  study,  we  consider  the  application  of 
vortex  merger  specifically  to  the  formation  of  hurricane  spiral  bands. 


Figure  3.24:  The  wave-activity  and  PV  flux  at  4  hours  intervals  for  experiment  two.  Times 
correspond  to  0  hrs  (upper  left),  4  hrs  (upper  right),  8  hrs  (lower  left),  12  hrs  (lower  right). 
Contours  are  from  5.0xl07  to  4.5  x  108  by  5.0xl07  and  are  scaled  by  10~7. 
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When  a  cyclonic  vortex  moves  sufficiently  close  to  another  cyclonic  vortex  of  equal 
intensity,  the  two  vortices  may  merge  or  coalesce  to  form  a  larger  vortex.  Extending  from 
the  newly  formed  vortex  will  be  spiral  arms  of  relatively  high  vorticity.  This  process  has 
been  shown  in  laboratory  experiments  by  Griffiths  and  Hopfinger  (1987).  If  the  vortices 
are  not  sufficiently  close  to  one  another,  they  may  simply  rotate  around  one  another 
without  merging.  This  is  sometimes  referred  to  as  the  Fujiwhara  effect  (Fujiwhara,  1923). 
Of  greater  importance  to  hurricane  spiral  bands  is  the  merger  of  vortices  having  different 
intensities.  We  consider  the  following  five  analogies  in  this  section,  all  of  which  produce 
spiral  bands  in  the  PV  field.  The  first  three  simulations  are  performed  using  a  transform 
grid  with  2562  points  and  keeping  108  waves.  Although  quadratically  nonlinear  terms  are 
not  alias  free,  results  indicate  that  for  relatively  short  model  integrations  (~1  day),  the 
increased  resolution  gained  outweighs  the  problem  of  nonlinear  computational  instability. 

The  first  simulation  involves  the  interaction  of  an  intense  vortex  with  a  weaker  vortex 
of  similar  size.  If  we  let  the  PV  associated  with  the  weak  vortex  represent  a  region  of  less 
organized  convection,  then  this  case  is  analogous  to  a  tropical  cylcone  which  encounters 
an  area  of  convection  on  the  scale  of  a  tropical  cloud  cluster.  We  present  this  simulation  in 
Fig.  3.25.  The  leftmost  vortex  is  initialized  with  a  maximum  vorticity  which  is  three  times 
that  of  the  rightmost.  As  the  vortices  evolve,  the  differential  rotation  associated  with  the 
more  intense  circulation  rapidly  acts  to  deform  the  weaker  vortex,  whereas  the  circulation 
associated  with  the  weaker  vortex  has  only  a  slight  impact  on  tue  more  intense  vortex.  By 
24  hours,  the  weaker  vortex  has  been  stretched  and  advected  towards  the  vortex  center 
yielding  a  single  spiral  band  of  PV. 

In  our  second  simulation,  we  consider  a  large,  strong  vortex  which  encounters  a 
smaller,  weaker  vortex.  This  case  is  analogous  to  a  tropical  cyclone  which  encounters 
a  relatively  small  region  of  convection.  Figure  3.26  shows  the  evolution  of  the  PV  field  for 
this  simulation.  The  intense  circulation  of  the  large  vortex  rapidly  deforms  the  smaller 
vortex.  By  24  hours,  the  PV  of  the  smaller  vortex  has  wrapped  around  the  larger  vor¬ 
tex  again  forming  a  thin  spiral  band.  We  should  therefore  also  expect  smaller  regions  of 
convection  in  nature  to  form  thinner  spiral  bands. 


two  initially  equally  sized  vortices  having  different  intensities.  The  leftmost  vortex  was 
initialized  with  a  maximum  relative  vorticity  which  was  three  times  that  of  the  rightmost. 
Units  are  10-5  s-1  with  a  contour  interval  of  10-4  s-1. 


Figure  3.26:  The  normalized  PV  field  at  0,  8,  16,  and  24  hours  showing  the  merger  of  a 
large,  strong  vortex  with  a  small,  weak  vortex.  The  large  vortex  was  initialized  with  a 
maximum  relative  vorticity  which  was  three  times  that  of  the  the  small  vortex.  Units  are 
10~5  s'1  with  a  contour  interval  of  10~4  s-1. 
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In  our  third  simulation,  we  consider  a  strong,  small  vortex  which  encounters  a  larger, 
weaker  vortex.  This  case  is  analogous  to  a  hurricane  which  forms  near  a  monsoon  trough 
which  is  rather  amorphous  in  appearance.  Figure  3.27  shows  the  evolution  of  the  PV 
field  for  this  simulation.  The  initial  vorticity  of  the  small  vortex  has  a  maximum  value 
which  is  three  times  that  of  the  larger  vortex.  As  time  evolves,  the  wind  field  associated 
with  the  small  intense  vortex  stretches  the  PV  field  of  the  large  vortex.  By  24  hours  the 
PV,  which  was  initially  associated  with  the  large  vortex,  completely  surrounds  the  smaller 
vortex.  In  addition,  the  remainder  of  the  larger  vortex  appears  as  a  wide  spiral  band 
extending  from  the  smaller  vortex.  The  large  region  of  PV  therefore  allows  the  smaller 
intense  vortex  to  grow  in  size.  This  may  help  explain  the  large  variation  in  size  of  tropical 
cyclones  (Merrill  ,  1984).  A  theory  has  recently  been  put  forth  by  G.  Holland  (personal 
communication)  which  suggests  that  if  a  tropical  cyclone  forms  near  the  monsoon  trough, 
it  would  be  provided  with  a  relatively  large  pool  of  PV  from  which  to  grow;  however,  if  a 
tropical  cyclone  forms  away  from  such  a  PV  source,  it  would  be  expected  to  remain  much 
smaller  in  size.  This  theory  will  also  be  examined  in  the  fourth  simulation. 

For  the  fourth  simulation,  we  consider  the  effect  of  placing  a  vortex  near  a  strip  of 
relatively  high  PV.  This  experiment  is  performed  using  1282  points  on  the  transform  grid 
while  keeping  42  waves,  and  is  intended  to  simulate  the  effect  of  tropical  cyclone  forming 
near  the  edge  of  a  monsoon  trough  or  ITCZ  which  is  zonal  in  appearance.  Figure  3.28 
shows  the  PV  field  for  this  simulation  at  0,  16,  32,  and  48  hours.  The  cyclonic  circulation 
associated  with  the  vortex  acts  to  advect  relatively  high  PV  air  from  the  strip  northward 
on  the  right  side  of  the  vortex  and  advect  relatively  low  PV  air  southward  on  the  left  side 
of  the  vortex.  As  time  evolves  the  vortex  is  made  larger  by  drawing  upon  the  PV  of  the 
strip.  By  48  hours  PV  from  the  strip  completely  surrounds  the  original  vortex  and  a  large 
spiral  band  has  formed.  This  result  again  supports  Holland’s  idea  for  the  development  of 
large  tropical  cyclones. 

In  the  previous  four  simulations  we  considered  the  formation  of  inner  bands  from 
initial  states  which  had  two  pre-existing  vortices  or  regions  of  PV.  For  the  fifth  and  final 
simulation,  we  consider  a  single  axisymmetric  vortex  which  encounters  an  area  of  intense 


Figure  3.27:  The  normalized  PV  field  at  0,  8,  16,  and  24  hours  showing  the  merger  of 
a  large  weak  vortex  with  a  small  strong  vortex.  The  small  vortex  was  initialized  with  a 
maximum  relative  vorticity  which  was  three  times  that  of  the  large  vortex.  Units  are  10~5 
s"1  with  a  contour  interval  of  10-4  s-1. 


Figure  3.28:  The  normalized  PV  field  at  0,  16,  32,  and  48  hours  showing  the  interaction  of 
a  vortex  with  a  strip  of  PV.  The  vortex  was  initialized  with  a  maximum  relative  vorticity 
of  4.0 x  10-4  s-1.  The  strip  was  initialized  with  a  maximum  relative  vorticity  of  6.25  x  10~5 
s-1.  Units  are  10~5  s_1  with  a  contour  interval  of  10-4  s_1. 
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convection.  The  goal  is  to  more  closely  simulate  the  actual  banding  process  by  convective 
asymmetries  in  nature. 

The  generation  of  PV  by  convective  heating  in  the  model  is  represented  as  a  mass  sink 
in  the  continuity  equation.  This  is  analogous  to  an  isentropic  (layers  of  constant  potential 
temperature)  model  in  which  heating  requires  mass  to  leave  a  layer  because  the  heated 
air’s  potential  temperature  changes,  thus  simulating  a  mass  transport  to  another  layer. 
Mathematically,  we  can  see  that  a  mass  sink  is  a  source  of  PV  by  examining  (3.5).  Since 
P  is  positive  for  cyclonic  circulations  in  the  northern  hemisphere,  a  mass  sink  (negative  S) 
will  produce  dP/dt  >  0  and  thus  be  a  source  of  PV.  The  conversion  between  a  mass  sink  in 
a  shallow  water  model  and  actual  heating,  however,  is  not  a  precise  one.  We  can  estimate 
an  appropriate  heating  rate  by  comparing  the  mass  source/sink  terms  of  the  continuity 
equations  for  both  a  shallow  water  system  and  a  stratified  fluid  in  isentropic  coordinates. 
That  is,  we  assume  pw  «  cr(),  where  a  —  —(l/g)(dp/d8),  w  is  the  vertical  velocity,  and  p 
is  the  density.  If  we  choose  149  kg/m2K  as  typical  value  of  a  for  a  tropical  atmosphere 
(Schubert  et  al.,  1992)  we  find  that  a  mass  sink,  w,  of  0.017  ms-1  corresponds  roughly  to 
a  heating  rate  of  10  K/day  (for  p  =1  kg/m3). 

Because  heating  is  added  to  the  model,  the  basic  state  depth  of  the  fluid  ( H )  is 
increased  to  2000  m,  which  allows  the  model  to  evolve  without  encountering  the  problems 
associated  with  massless  layers.  By  increasing  the  basic  state  depth  of  the  fluid,  however, 
we  also  increase  the  pure  gravity-wave  speed  and  therefore  the  Rossby  radius.  Model 
results  indicate  that  these  changes  have  no  adverse  effects  on  the  simulations. 

Our  simulation  consists  of  an  axisymmetric  vortex  which  is  located  north  of  a  large 
convection  area.  This  area  of  convection  is  an  ellipse  with  tin  aspect  ratio  of  two  whose 
center  is  initially  located  400  km  from  the  center  of  the  vortex.  The  elliptical  region  has 
a  maximum  semi-major  axis  length  of  300  km  and  the  x  and  y  dependence  of  the  heating 
uses  the  Melander  profile  with  a  steepness  parameter  of  2.56,  i£,=100  km  and  Ro=300 
km.  The  heating  function  is  therefore  similar  in  shape  (but  not  size)  to  the  initial  vorticity 
pattern  for  experiment  two  in  section  3.4.  The  maximum  value  of  the  heating  function  is 
0.08  ms-1  which  equates  to  a  maximum  heating  rate  of  ~40-50  K/day. 
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Figure  3.29:  The  time  dependent  part  of  the  heating  function  Q(x,y,t). 

To  prevent  the  excitation  of  transient  gravity  waves  and  ensure  evolution  along  the 
diabatic  slow  manifold  (Schubert  et  a 1980),  the  heating  is  turned  on  slowly  using  a 
function  of  the  form 

Q(x,  y,  t)  =  Q(x,  y)[l  -  exp(-oi2)],  (3.42) 

where  is  the  constant  a  is  3.5xl0~9s-2  and  Q(x,  y)  is  the  elliptical  shaped  heating  pattern 
described  above.  The  constant  a  was  chosen  such  that  the  heating  reaches  its  peak  at 
exactly  12  hrs.  The  model  is  integrated  using  1282  points  on  the  transform  grid  while 
keeping  42  waves.  Fig.  3.29  shows  how  (3.42)  varies  with  time. 

Results  from  this  simulation  are  shown  in  Figs.  3.30-3.32.  Initially  the  vortex  was 
located  200  km  from  the  top  of  the  domain  and  centered  in  x.  The  heat  source  is  stationary 
and  located  200  km  from  the  bottom  of  the  domain  and  also  centered  in  x.  By  24  hrs,  the 
presence  of  the  heating  anomaly  is  felt  on  the  entire  right  edge  of  the  vortex.  In  addition, 
the  vortex  has  shifted  towards  the  lower  pressure  created  by  the  heating,  moving  almost 
200  km.  The  weak  cyclonic  circulation  created  by  the  heating  anomaly  has  caused  the 
vortex  also  to  move  slightly  westward.  At  36  hrs,  the  PV  anomaly  produced  by  the  heating 
has  stretched  and  shifted  to  encompass  the  entire  northern  half  of  the  vortex.  By  48  hrs, 
the  PV  field  from  the  vortex  and  that  due  to  the  heating  have  merged.  The  result  is  a 
large  tail  of  PV  extenting  in  a  spiral  fashion  from  the  southern  edge  of  the  vortex.  The 
height  field  associated  with  the  vortex  shows  slight  asymmetries  associated  with  the  spiral 


78 


Figure  3.30:  The  wind  and  mass  fields  for  the  asymmetrically  heated  vortex  at  24  hrs.  The 
height  (upper  right)  is  measured  in  decameters  and  has  a  contour  interval  of  2  decameters. 
The  maximum  wind  vector,  defined  as  the  distance  between  two  consecutive  tick  marks, 
is  50  ms-1.  The  normalized  PV  (lower  left)  and  absolute  vorticity  (lower  right)  are  in 
units  of  10-5  s_1  and  have  contour  intervals  of  2.0xl0-4  s-1.  Shown  is  the  inner  800x 
800  km  domain.  Tick  marks  are  every  50  km. 
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Figure  3.31:  The  same  as  Fig.  3.30  but  for  36  hrs. 
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Figure  3.32:  The  same  as  Fig.  3.30  but  for  48  hrs. 
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band  suggesting  the  model  is  evolving  on  the  diabatic  slow  manifold.  The  band  is  also 
found  to  correspond  to  the  location  of  a  secondary  tangential  wind  maximum.  Figure 
3.33  shows  the  tangential  wind  speed  for  four  radial  arms  which  emanate  from  the  vortex 
center  at  90°  angles  starting  from  the  positive  i  axis.  The  upper  left  panel  shows  clearly 
that  the  band  is  associated  with  an  increase  in  the  tangential  wind.  These  secondary  wind 
maxima  are  often  found  in  nature  (Willoughby  et  a/.,  1982). 

These  results  seem  to  indicate  that  banding  process  can  be  explained  rather  well  using 
PV  arguments.  Other  related  results  (not  shown)  indicate  that  asymmetric  heating  may 
also  be  responsible  for  the  formation  of  concentric  eyewalls  which  are  observed  occasionally 
in  smaller  intense  tropical  cyclones  (Willoughby  et  ad ,  1982).  During  a  model  simulation  in 
which  the  heating  was  turned  off  slowly,  thus  simulating  a  vortex  which  comes  in  contact 
with  an  area  of  intense  convection  but  slowly  propagates  away,  the  asymmetries  in  the  PV 
field  simply  wrapped  around  the  existing  vortex,  creating  an  annular  region  of  PV.  The 
fact  that  concentric  eyewalls  are  only  observed  in  the  more  intense  hurricanes  supports 
this  idea.  When  coming  in  contact  with  an  area  of  intense  convection,  the  more  intense 
hurricanes  would  allow  the  PV  anomalies  to  simply  wrap  around  the  vortex  before  having 
time  to  interact  with  it.  Observational  studies  would  need  to  be  done,  however,  to  support 
this  hypothesis. 
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Figure  3.33:  The  tangential  wind  speed  along  four  radial  arms  extending  from  the  center 
of  the  vortex  at  48  hrs  (see  Fig.  3.32).  The  radial  arms  sure  at  0°  (lower  left),  90°  (lower 
right),  180°  (upper  left),  and  270°  (upper  right)  as  measured  from  the  positive  x  axis. 
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OUTER  BAND  FORMATION  THROUGH  ITCZ  BREAKDOWN 

In  this  chapter  we  examine  the  formation  of  outer  bands  in  tropical  cyclones  result¬ 
ing  from  the  breakdown  of  the  intertropical  convergence  zone  (ITCZ).  We  begin  by  first 
defining  outer  bands  and  presenting  observational  evidence  for  their  existence.  In  section 
4.2  we  present  a  brief  review  of  the  applicable  general  stability  theories  which  help  us  to 
understand  the  necessary  conditions  for  barotropic  instability.  These  theories  provide  a 
basis  for  understanding  some  of  the  fundamental  dynamics  of  ITCZ  breakdown.  In  the 
final  section,  we  present  both  numerical  analysis  and  model  simulations  of  ITCZ  break¬ 
down  and  outer  band  formation.  The  numerical  analysis  makes  use  of  Unear  normal  mode 
theory  which  (although  it  is  sfightly  less  general  than  other  stabihty  theories)  allows  us  to 
predict  unstable  patterns  and  growth  rates.  The  normal  mode  stabihty  analysis  is  then 
foUowed  by  five  numerical  experiments. 

4.1  Observational  Evidence  for  Outer  Bands 

Before  presenting  evidence  for  the  existence  of  outer  bands,  it  is  first  necessary  to 
define  exactly  what  is  meant  by  outer  bands.  Outer  hurricane  bands  are  defined  as  bands 
of  convection  which  are  typically  located  farther  than  500  km  from  the  center  of  the  vortex 
and  owe  their  existence  to  the  breakdown  of  the  ITCZ.  As  an  example,  we  present  GOES 
IR  satelhte  images  of  the  tropical  North  Pacific  during  late  July  1988  in  Figs.  4. 1-4.2. 

On  July  26,  the  ITCZ  appears  as  a  relatively  uniform  area  of  convection  which  varies 
between  5°  and  9°  latitude  in  width  and  nearly  40°  latitude  in  length.  By  August  3rd, 
the  ITCZ  has  broken  down  and  five  tropical  disturbances  remain.  Between  some  of  the 
disturbances  there  exist  noticeable  filaments  of  convection.  These  are  the  outer  bands. 
By  August  12th,  the  ITCZ  has  again  begun  to  reform,  and  the  process  begins  to  repeat 


Figure  4.1:  GOES  IR  images  at  1646  UTC  on  26  July  (top)  and  3  August  (bottom),  1988. 
Outer  bands  are  especially  noticeable  extending  from  the  second  vortex  from  the  left  in 
the  bottom  panel. 


. mis . i . . . . 

Figure  4.2:  Same  as  Fig.  4.1  but  for  12  August  1988. 

itself.  The  basic  dynamics  of  the  ITCZ  breakdown  process  and  outer  band  formation  can 
be  explained  through  barotropic  instability  arguments.  This  is  not  to  imply  that  moist 
physical  processes  and  baroclinicity  do  not  play  a  role.  Certainly,  they  must.  We  are 
simply  trying  to  isolate  the  fundamental  aspects  of  the  dynamics. 

A  relationship  between  the  ITCZ  and  hurricane  bands  was  first  noted  by  Fletcher 
(1945).  Fletcher  hypothesized  that  spiral  bands  were  formed  by  individual  convergence 
lines  from  the  ITCZ  being  “coiled”  into  the  center  of  the  tropical  storms.  Although  his 
theory  is  flawed,  the  connection  he  drew  between  the  ITCZ  and  hurricane  bands  may  be 
valid. 

In  the  following  two  sections,  we  will  show  that  the  breakdown  of  a  strip  of  relatively 
high  potential  vorticity  (PV)  air  takes  on  many  of  the  same  attributes  of  the  observed  ITCZ 
breakdown,  if  regions  of  higher  PV  are  interpreted  as  regions  of  more  intense  convection. 
The  breakdown  of  a  PV  strip  is  characterized  by  the  “pooling”  of  PV  in  certain  locations 
which  causes  other  regions  of  PV  to  become  extremely  thin  and  eventually  detach  ( Pratt 
and  Pedlosky.  1991).  The  mechanism  for  this  breakdown  was  best  described  by  Hoskins 


86 


Table  4.1:  Summary  of  linear  and  nonlinear  stability  theory  for  various  types  of  dynamical 
systems.  Question  marks  imply  no  known  theory  exists. 


Dynamics 

Linear 

Nonlinear 

Nondivergent 

Rayleigh  1880 

Taylor  1915 

Fj0rtoft  1950 

Amol’d  1966 

Quasi-geostrophic 

Chamey-Stem  1962 

Shepherd  1988 

Semi-geostrophic 

Eliassen  (/-plane)  1983 

Magnusdottir  (,0-plane)  1990 

Magnusdottir  (sphere)  1991 

? 

Primitive  Equation 

Ripa  (shallow  water)  1983 

Ripa  (multiple  finite  layers)  1991 

? 

efc  aJ.  (1985)  as  the  cooperative  interaction  of  counter  propagating  Rossby  waves,  or  more 
generally,  PV  waves.  We  begin  our  discussion  by  first  familiarizing  ourselves  with  the 
basic  principles  of  barotropic  stability  theory. 

4.2  A  Review  of  Stability  Theory 

In  this  section,  we  review  some  of  the  fundamental  barotropic  stability  theorems. 
These  theorems  offer  insight  to  the  understanding  of  ITCZ  breakdown  by  providing  the 
conditions  necessary  for  barotropic  instability.  The  present  state  of  stability  theory  for 
various  dynamical  systems  is  summarized  in  Table  4.1.  Currently,  nonlinear  stability 
theorems  are  restricted  to  highly  balanced  flows.  Here  we  review  the  existing  stability 
theorems  for  both  the  nondivergent  barotropic  equations  and  the  shallow  water  (primitive) 
equations.  We  begin  by  reviewing  Rayleigh’s  famous  stability  theorem  for  the  linear 
nondivergent  barotropic  system.  This  discussion  is  followed  by  stability  theorems  for  the 
nonlinear  nondivergent  system  as  derived  by  Amol’d  (1965,66)  and  the  linearized  shallow 
water  system  as  derived  by  Ripa  (1983). 

4.2.1  Rayleigh’s  Theorem 

The  most  elegant  method  of  deriving  Rayleigh’s  stability  theorem  was  first  presented 
by  Taylor  (1915).  This  method  involves  the  use  of  the  variable  t]  which  represents  a  particle 
displacement  in  the  y-direction.  We  begin  the  derivation  by  considering  the  nondivergent 
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vorticity  equation  in  Cartesian  coordinates  linearized  about  a  basic  state  which  has  only 
a  zonal  component.  The  absolute  vorticity  (  here  defined  as  0  equation  and  particle 


displacement  equations,  take  the  form 


£*•!- 


3-' 

where  D/Dt  =  d/dt+udjdx.  The  barred  quantities  represent  zonal  means  whereas  primed 
quantities  denote  small  amplitude  perturbations  about  the  zonal  mean.  By  combining 


(4.1)  and  (4.2)  we  obtain 


l  M)- 


If  the  quantity  inside  the  parentheses  in  (4.3)  vanishes  at  t  —  0,  then  it  vanishes  for  all 
time.  Multiplying  this  quantity  by  —v'  and  rearranging  terms  gives  us  the  relationship 

££  -  v'f  =  0,  (4.4) 


where 


A  s  -4„2|£ 


represents  the  wave  activity  for  the  nondivergent  barotropic  system.  If  we  expand  C  and 
use  the  continuity  equation,  (4.4)  can  be  written  in  final  form  as 

dA  ,  d[uA  +  5(u'2  -  u'2)]  _  0(uV)  n 

ir+ - ^ - +_*r  0  (46) 

We  note  here  the  similarity  between  the  linear  wave-activity  equation  (4.6)  and  the  non¬ 
linear  shallow  water  form  (3.39). 

Rayleigh’s  stability  theorem  can  now  be  obtained  by  integrating  (4.6)  over  the  entire 
domain  (assuming  a  doubly  periodic  domain)  to  obtain 


7tIIAixdy  =  0- 


From  (4.5)  and  (4.7),  we  see  that  if  the  mean  vorticity  field  is  monotonic,  t/  can  grow  only 
locally.  If  77  is  to  grow  in  an  overall  sense,  the  mean  vorticity  gradient  must  reverse  sign. 
This  is  Rayleigh’s  classic  necessary  condition  for  barotropic  instability. 
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4.2.2  Arnol’d’s  Theorem 

In  the  previous  subsection  we  considered  the  linear  nondivergent  barotropic  governing 
equations.  We  now  wish  to  generalize  Rayleigh’s  theorem  to  finite  amplitude  disturbances. 
This  was  first  accomplished  by  Amol’d  (1965, 1966).  In  this  subsection,  we  derive  Amol’d’s 
theorem  for  the  /-plane  in  Cartesian  coordinates.  The  governing  equations  for  this  case 
can  be  written  in  the  form 

-  fiy  - !/«)]  +  -  f(y  -  vc)) + p/p)  +  ~  /(»  -  yc))]  =  o  (4.8) 

dv  d(uv)  dk(vv)  1  dp  .  . 

9f+-fei+'V'+/“+;f  ,4-9) 


du  dv  _ 
dx  +  dy 


(4.10) 


where  yc  is  a  constant  to  be  defined  later.  We  regard  (4.8)-(4.10)  as  a  closed  system  in  u,  v 
and  p,  where  the  density,  p,  is  constant.  A  consequence  of  (4.8)-(4.10)  is  the  conservation 
of  the  absolute  vorticity  £.  Since  (  is  conserved,  any  function  C(C)  is  also  conserved. 
Thus,  the  Casimir  function  C(Q  satisfies 


dc  d(uC)  d(vC) 

dt  dx  dy 


(4.11) 


Adding  (4.8)  and  (4.11),  we  obtain 


d[«  ~  /(y  -  yc)  +  C]  a[u(u  -  f(y  -  ye)  +C)+  p/p)  a[u(u  -  f{y  -  yc)  +  C)]  _ 


(4.12) 

Following  McIntyre  and  Shepherd  (1987)  and  Haynes  (1988),  we  now  wish  to  choose  C(() 
such  that  the  disturbance  part  of  [u  —  f(y  -  yc)  +  C]  can  be  written  as  a  divergence, 
plus  a  term  which  is  explicitly  second  order  in  wave  amplitude.  The  disturbance  part 
of  [u  -  /(y  -  yc)  +  C]  is  defined  as  (u  -  /(y  -  yc)  +  C]e  =  ue  +  C«)  -  C(Co)-  The 
subscript  e  is  now  used  to  represent  the  disturbance  part  of  a  field  while  a  subscript  0 
represents  its  basic  state.  This  change  in  notation  is  to  distinguish  the  finite  amplitude 
disturbances  discussed  in  this  subsection  from  the  small  amplitude  approximation  of  the 


previous  subsection.  After  some  manipulation,  the  above  expression  can  be  written 


(u  -  /(y  -  yc)  +  C]e  =  A  + 


flJC'KoKI  <?[C'«oKJ 


dC«o) 


(4.13) 


where 


A  =  C(()  —  C(Co)  -  CeC (Co) 

is  the  wave-activity  of  the  system  and  primed  quantities  now  refer  to  derivatives  with 
respect  to  the  dependent  variable.  Through  use  of  a  Taylor  series  expansion,  it  is  easy  to 
verify  that  A  is  a  second  order  quantity.  In  addition,  the  second  term  on  the  right  hand 
side  of  (4.13)  is  the  divergence  of  a  flux.  We  now  wish  to  choose  C(C)  in  such  a  way  that 


the  last  term  in  (4.13)  vanishes.  The  choice 


C(0  =  j^[yo(0  ~  Vc]dC, 


(4.14) 


where  (c  =  (o(t/c)  and  yo(0  is  the  inverse  of  Co(y)  (i.e.  yo(Co(y))  =  y),  makes  the  last  term 
in  the  square  brackets  in  (4.13)  vanish,  as  can  easily  be  verified. 

Noting  that  the  basic  state  is  a  steady  solution  of  the  governing  equations  and  then 
substituting  (4.13)  into  (4.12)  we  obtain 


~di  +dx  [°^U  ~  +  ^  +P/P  + 

+!  [»(<•  -  /(»  -  yc)  +  C)  +  -CKo)—-]  =  0. 


(4.15) 


Using  the  disturbance  momentum  equations 

fae  ,  .  d(%uZ  +  uoue  +  %v*+pe/p) 

-ar-Cv'  + - % - =0’ 

aT  +  <u  "  +  ^ ^ - =  o 

in  (4.15)  and  defining  yc  such  that  uo(yc)  =  0,  we  obtain 

dA  d  .  .  t,2  2m  3 .  .  .  n 

~dt  +  dx  uA  +  ^Ue  ~  +  dy^VeA  +  UeVe ^  =  °’ 

where  the  wave-activity  can  also  be  written  as 


(4.16) 


(4.17) 


(4.18) 


A  =  -  [  [yo«o  +  0  -  yo«o)]d<. 

Jo 


(4.19) 


To  obtain  the  nonlinear  stability  condition  we  now  integrate  (4.18)  over  the  entire 


domain  (again  assuming  a  doubly  periodic  domain)  to  obtain 


±JjAdxdy  =  0. 


(4.20) 
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Although  this  looks  identical  to  (4.7),  we  must  keep  in  mind  that  the  A  in  (4.20)  is  defined 
by  (4.19)  while  the  A  in  (4.7)  is  defined  by  (4.5).  The  results  are  consistent,  however,  since 
(4.19)  reduces  to  (4.5)  in  the  small  amplitude  limit.  To  see  this,  we  consider  the  integrand 
in  (4.19),  i.e. 

[yo(Co  +  0  -  yo(Co)]  (4.21) 

A  Taylor  series  expansion  about  Co  of  (4.21)  yields  (dyo(Co)/dC)C-  Using  this  result  to 
evaluate  the  integral  in  (4.19)  yields 

A  =  -<'2^r^),  (4.22) 

where  (e  =  ('  in  the  small  amplitude  limit.  Using  the  relationship  given  by  (4.3)  to 
eliminate  ('  and  noting  that  dyo/d(  =  (d(o/dy)~l  =  (d£/dr)~l  in  the  small  amplitude 
limit,  we  arrive  at  the  desired  relationship.  That  is, 

clim  j-  jf  [yo(Co  +  0  ~  yo(Co)]«*C  j  =  ~W fy-  (4-23) 

Thus,  (4.19)  is  simply  a  finite  amplitude  generalization  of  the  wave-activity  obtained 
through  Rayleigh’s  theorem. 

We  now  note  that  since  yo  Is  a  monotone  function  of  C,  the  integrand  in  (4.19), 
is  bounded  by  the  lines  — CMyo/dC|max  aud  — C|dyo/dC|min-  These  lines  represent  the 
maximum  and  minimum  slopes  respectively  of  the  integrand  in  (4.19)  multiplied  by  the 
variable  of  integration,  (.  In  addition,  the  monotonic  nature  of  yo  implies  A  <  0  for  all 
values  of  Ce-  This  allows  us  to  write 

C«2Myo/dC|min  <  2|A(C0,Ce)|  <  CVyo/dCImax.  (4.24) 

Together,  (4.20)  and  (4.24)  imply  that 

!<fyo/dCUin  Jf  C?(x,y,t)dxdy  <  JJ 2|A(C0,Ce,t)|  dxdy  = 

yj2lA(Co,Ce>0)l  dxdy  <  |dy0/d<|max  J f  Ce(*-y.O )dxdy, 
which  can  also  be  written 


A 
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This  is  the  form  of  Arnol’d’s  (1965,  1966)  result  derived  by  McIntyre  and  Shepherd  (1987). 
The  inequality  (4.25)  bounds  the  disturbance  enstrophy  at  time  t  in  terms  of  the  initial 
disturbance  enstrophy  and  the  radial  gradient  of  the  basic  state  absolute  vorticity.  It 
rules  out  the  possibility  of  instability  for  basic  state  flows  with  d^o/dy  >  0  (or  d^o/dy  < 
0)  everywhere.  Thus,  Amol’d’s  theorem  provides  a  sufficient  condition  for  stability,  or 
conversely,  a  necessary  condition  for  instability. 

4.2.3  Ripa’s  Theorem 


Until  now,  we  have  discussed  stability  with  regard  only  to  nondivergent  flow.  We  now 
wish  to  relax  this  assumption  of  balance  but  retain  the  assumption  of  linearity.  The  case 
we  consider  is  the  linear  shallow  water  equations  on  the  /-plane.  This  theorem  was  first 
derived  by  Ripa  (1983)  on  the  /3-plane  and  sphere  for  a  single  layer  of  fluid.  Ripa  (1991) 
later  generalized  his  theorem  to  include  multiple  discrete  incompressible  layers  (but  not 
continuously  stratified). 

The  shallow  water  equations  linearized  about  a  basic  state  which  is  in  gradient  balance 
can  be  written  in  cylindrical  coordinates  as 


where  P  = 
we  need  to 


dv!  -  ,  _  d(uu'  +  gh')  A 

q,  '  a  U, 

at  ox 

g  +  fa  +  PM  +  W+m,  o, 

dh!_  d(hu'  +  uh')  d(i/h) 

dt  dx  dy  ’ 


(4.26) 

(4.27) 

(4.28) 


C/h,  C  being  the  zonally  averaged  absolute  vorticity.  To  derive  Ripa’s  theorem 
combine  (4.26)-(4.28)  into  equations  for 


Ef  =  \  ]h(v!2  +  v'2)  +  2 uu'ti  +  gh'2}  , 

(4.29) 

M'  =  u'h', 

(4.30) 

and 


(4.31) 
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To  derive  the  equation  for  E 1  we  form  (hu'  +  h'ii)  •  (4.26)  ■+■  (hv1)  •  (4.27)  4-  (gh1  +  uu')  ■  (4.28) 
to  obtain 


dE1  r2_  fp/  d[(hu'  4-  uh')(gh!  +  uu'))  d[hv'(uu'  +  gh')] 
&t  UV  dx  dy 

To  derive  the  equation  for  M'  we  form  h!  ■  (4.26)  +  u'  •  (4.28)  to  obtain 

dM'  p uipi  +  OfiM*  4-  lh(u*  -  v'2)  +  +  dfttVh] 

dt  dx  dy 


(4.32) 


(4.33) 


Combining  the  integrated  forms  of  (4.32)  minus  uo  times  (4.33)  we  obtain 


d(E'  -  uoM') 
dt 


+  h2v'P'{a  -  uo)j 


dxdy  as  0, 


(4.34) 


where  uo  is  an  arbitrary  constant.  The  equation  for  P'  is  obtained  by  forming  the  vorticity 
equation  from  (4.26)-(4.27)  and  eliminating  the  divergence  using  (4.28).  This  equation 
takes  the  form 


DP1  dP  ,  n 

dT  +  ^"  =0’ 


(4.35) 


where  D/Dt  =  dfdt  +  ud/dx.  Defining  the  meridional  particle  displacement,  r/,  by 
Dy/Dt  —  v1  and  assuming  both  P*  and  rf  vanish  at  some  finite  time,  we  obtain 


^  dP  A 
P'  +  —r)  =  0. 

dy 


(4.36) 


Multiplication  of  (4.36)  by  v'  yields 


(4.37) 


Finally,  using  (4.37)  in  (4.34)  for  v'P1  we  obtain 


d_ 

dt 


jj  ^E/-  uo M'  -  dxdy  =  °i 


(4.38) 


which  is  the  divergent  barotropic  generalization  of  the  nondivergent  bairotropic  result  (4.7). 
We  now  argue  that,  if  E'  —  wqM'  >  0  and  (u  —  uo )dP/dy  <  0,  the  constraint  (4.38)  does 
not  allow  rj2  to  grow  in  an  overall  sense.  We  now  note  that  the  three  rightmost  terms 
(4.29)  can  be  written  in  a  quadratic  form  for  h!  and  v' .  The  quadratic,  and  therefore  E1, 
will  be  positive  if  the  coefficients  of  the  quadratic  satisfy  h  >  0  and  u2  <  gh  (Strang,  1988, 
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p.  325).  Since  h  >  0  by  assumption,  then  E?  >  0  if  u2  <  gh.  By  a  similar  argument, 
E *  —  wo M1  >  0  if  [(u  —  tio)]2  <  gh.  We  can  now  state  Ripa’s  shallow  water  generalization 
of  the  theorems  of  Rayleigh  and  Fjortoft.  If  there  exists  any  value  of  uq  such  that 

JD 

(ti  —  uq)—  <  0  and  (ti  —  uo)2  <  gh  (4.39a, b) 

ay 

for  all  r,  then  the  flow  is  stable  to  infinitesimal  perturbations.  Ripa  has  also  discussed 
several  corollaries  of  (4.39),  one  of  which  is  obtained  by  choosing  uo  =  max[uj.  This 
results  in  the  following  weaker  sufficient  condition  for  stability.  If 

—  >0  and  max[ti]  <  min[ti  +  {gh)1/2]  (4.40a, b) 

dy 

for  all  r ,  then  the  flow  is  stable  to  infinitesimal  perturbations. 

To  recover  the  stability  results  for  the  nondivergent  barotropic  model  from  the  sta¬ 
bility  results  for  the  divergent  barotropic  model  we  consider  the  limit  gh  — »  oo,  in  which 
case  (4.39b)  is  satisfied  for  any  finite  tio*  Then,  there  is  no  difference  between  vorticity 
and  potential  vorticity,  and  a  choice  of  tio  such  that  u  —  «o  <  0  everywhere  leads  to 
dC,/dy  >  0  everywhere  as  sufficient  for  stability,  while  a  choice  of  uo  such  that  u  —  uo  >  0 
everywhere  leads  to  d^/dy  <  0  everywhere  as  sufficient  for  stability.  Thus,  a  necessary 
condition  for  instability  is  that  d£/dy  have  both  signs  (Rayleigh’s  theorem).  It  is  also  of 
interest  to  note  that,  if  d^fdy  =  0  at  y  =  y,  then  the  choice  uo  =  u(y)  leads  from  (4.39a) 
to  (u(y)  —  uo]d(/dy  <  0  somewhere  as  a  necessary  condition  for  instability  (Fjortoft’s 
theorem). 

4.3  Numerical  Analyses  and  Simulations 

In  the  previous  section,  we  saw  that  the  common  factor  to  all  stability  analyses  was 
the  need  for  a  reversal  in  the  PV  gradient.  They  did  not  however  indicate  to  what  extent 
these  disturbances  could  grow  nor  with  what  speed.  In  this  section,  we  attempt  to  shed 
some  light  on  this  problem  through  both  normal  mode  stability  an<*lysis  and  numerical 
simulations.  This  theory  is  then  applied  to  the  breakdown  of  the  ITCZ  and  the  formation 
of  outer  bands. 
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In  the  first  subsection  we  present  a  linear  normal  mode  stability  analysis  of  a  PV  strip. 
This  provides  an  estimate  of  both  unstable  patterns  and  growth  rates  In  subsection  4.3.2 
we  simulate  the  breakdown  of  the  ITCZ  and  the  formation  of  outer  bands  through  five 
mouel  integrations. 

4.3.1  Normal  Mode  Stability  Analysis 

In  the  previous  section  we  examined  various  general  stability  theories.  Although  these 
theories  provide  necessary  conditions  for  stability,  they  tell  us  nothing  about  the  rate  at 
which  these  instabilities  will  grow  or  the  favored  patterns  for  instability.  To  determine 
these  unknowns,  we  utilize  the  less  general  tool  of  normal  mode  stability  analysis.  This 
type  of  analysis  is  less  general  in  the  sense  that  it  assumes  wave  form  instabilities,  i.e. 
stabilities  which  are  proportional  to  The  normal  mode  stability  analysis  of  a 

PV  strip  was  first  presented  by  Rayleigh  (1945,  vol.  2,  pp.  384-398)  and  is  presented  here 
as  a  review.  Note  that  for  the  nondivergent,  barotropic  analysis,  there  is  no  difference 
between  PV  and  absolute  vorticity. 

We  begin  the  analysis  by  considering  the  linearized  nondivergent  barotropic  dynamics 
of  a  PV  field  with  the  following  basic  state  zonal  flow 

{-Coy o  ya  <  y  <  °° 

-Coy  -y0<y<y0  .  (4.41) 

C oy0  -oo  <  y  <  -y0 

The  basic  state  absolute  vorticity  fiek»  corresponding  to  (4.41)  is 

/  y0  <  y  <  oo 

/  +  Co  -y0<y<y0  (4.42) 

/  -oc  <  y  <  -ya 

Because  the  vorticity  is  constant  everywhere  except  at  the  interfaces,  the  equation  which 
the  perturbation  streamfunction  must  satisfy  is  simply 

d V  ay 

dx 2  5y2 


«»>-'- 1  =  < 


(4.43) 
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If  we  assume  a  normal  mode  solution  which  is  proportional  to  e*(kx  vt\  the  solution  to 
(4.43)  can  be  written  as 

'  yne-k(v-v.)  +  yj^y+vo)  y0<y<oo 
tp'(x,y,t)  =  ei('kx~vt^  <  #ne«v-v.)  +  *7^+">>  -y0  <  y  <  y„  -  (4.44) 

(  -oo  <  y  <  -Vo 

It  is  easily  confirmed  by  inspection  that  ip*(x,  y,  t)  is  continuous. 

Now  let  us  assume  that  the  northern  and  southern  interfaces  of  the  PV  strip  are 
perturbed  in  a  sinusoidal  fashion  by  an  amount  r/„  =  Tjn^kx~l/t)  and  =  rjse,(kx~t/t), 
where  yn  and  tjs  represent  the  distances  the  northern  and  southern  interfaces  are  displaced 
from  their  original  positions  respectively.  Since  PV  is  materially  conserved,  tj„  and  ya 
also  represent  particle  displacements.  The  equations  governing  the  displacement  of  the 
boundaries  can  be  written 


=  t/  =  ^ 

dt  dx  8x 


dru  -&h  ,  _  <W_ 

dt  U  dx V  dx' 


(4.45) 


Our  goal  is  to  use  (4.44)  in  (4.45)  to  determine  under  what  conditions  the  particle  dis¬ 
placement  will  grow  unbounded.  Before  this  can  be  done,  however,  we  must  first  relate 
the  complex  streamfunction  amplitudes  (¥„  and  ¥,)  to  the  complex  particle  displacement 
amplitudes  ( »)„  and  ya).  We  do  this  by  requiring  the  zonal  component  of  the  wind  (u-f  u') 
to  be  continuous  at  the  boundaries  y  =  y0  +  rjn  and  y  =  —y0  +  ijt,  where  the  the  zonal 
wind  is  obtained  from  u1  —  —chp'/dy.  In  doing  this,  we  obtain 


=  $<<,»)„  and  k*s  =  -£< oVs- 


(4.46) 


Substituting  (4.46)  into  (4.44)  and  using  the  result  in  (4.45)  gives 

/  5<o(l-2fcy0)  -%(0e~2ky°  \  /  Vn  \  _  j 

\  %(0e-2kv*  -^C„(l-2 ky0)  )  \  fjs  )  ' 


(4.47) 


Equation  (4.47)  is  a  standard  matrix  eigenvalue  problem  which  can  be  solved  in  the  usual 
manner.  The  solutions  for  the  eigenvalues/frequencies  are 

i/fc  =  ±£<0[(l-2*y0)2-e-4fc*]5 


(4.48) 
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The  solution  for  v  will  be  purely  imaginary  if  ky0  lies  between  0  and  0.6392.  By  differen¬ 
tiating  (4.48)  with  respect  to  kya  and  setting  the  result  to  zero,  the  most  unstable  mode 
is  found  to  occur  when  kya  =  0.3984.  Since  k  =  2ir/Lx  where  Lx  is  the  wavelength  in  the 
x-direction,  we  obtain  the  following  relationship  between  the  half-width  of  the  disturbance 
and  the  wavelength  of  the  most  unstable  mode, 

2ny0  «  0.4 Lx,  (4.49) 

i.e.,  the  wavelength  of  the  most  unstable  mode  is  approximately  eight  times  the  width  of 
the  strip.  Thus,  in  the  absence  of  any  adverse  shear  or  variation  of  the  earth’s  vorticity, 
we  would  expect  the  wavelength  of  the  most  unstable  mode  to  be  constant  regardless  of 
the  intensity  of  vorticity  strip. 

The  above  analysis  can  be  made  slightly  more  general  by  including  the  effects  of 
adverse  shear,  that  is,  the  addition  of  a  uniform  anticyclonic  shear  layer  (Dritschel,  1989). 
We  include  this  effect  partly  because  of  generality  but  mainly  out  of  necessity.  The 
necessity  of  this  analysis  stems  from  the  models  periodic  boundary  conditions.  Having 
periodic  boundary  conditions  implies  no  mean  vorticity  field  can  exist  since  this  would 
imply  a  mean  circulation  around  the  periphery  of  the  model  domain.  A  mean  circulation, 
however,  requires  the  wind  to  be  moving  in  opposite  directions  at  opposing  boundaries. 
This  clearly  violates  the  models  requirement  of  periodicity.  As  will  be  seen  in  the  following 
section,  the  mean  vorticity  field  can  be  removed  by  subtracting  it  from  the  initial  vorticity 
field.  Subtracting  the  mean  positive  vorticity  is  equivalent  to  adding  a  uniform  region  of 
negative  vorticity.  We  denote  the  negative  vorticity  as  being  some  fractional  amount  of 
(o,  namely  —  AC,  where  A  can  be  any  number.  The  resulting  mean  wind  and  absolute 
vorticity  profiles  for  this  case  take  the  form 

{-< o(yo  -  Ay)  ya  <  y  <  oo 

-<<>y(  1  -  A)  -y0<y<y0  ,  (4.50) 

Co{yo  +  Ay)  -oo  <  y  <  -ya 
and 

(  f  -  CoA  ya  <  y  <  oo 


/  +  CoU  —  A)  -yQ  <y  <y0 
I  -  CoA  -oo  <  y  <  -y0 


(4.51) 
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Because  the  vorticity  is  constant  except  near  the  boundaries,  t/>'  must  still  satisfy  (4.44). 
In  addition,  (4.46)  remains  valid.  The  effect  of  adverse  shear  is  felt  only  in  the  definition 
of  u(y)  and  therefore  only  in  the  equations  for  rjn  and  rja.  Using  this  fact,  we  see  that  the 
new  matrix  eigenvalue  problem  takes  the  form 

/  jCol1  ~  2ky0(l  -  A)]  -\C0e~2kv“  \  /  4.  \  /  Vn  \ 

V  hCoe~2ky-  -Ko[l  ~  2fcy0(l  -  A)]  )  {  fjs  )  ~  V  \  Vs  )  '  ' 

The  corresponding  solution  to  (4.52)  is 

|*(A)  =  ±fa{[  1  -  2fcyo(l  -  A))2  -  e-4fc*>]5.  (4.53) 


Clearly  instabilities  cannot  exist  if  A  >  1  since  the  radicand  in  (4.53)  will  be  positive  for 
all  wavelengths.  This  will  result  in  v  being  strictly  real.  The  effect  of  adverse  shear  is  to 
decrease  both  the  wavelength  of  the  most  unstable  mode  and  its  growth  rate.  Figure  4.3 
shows  the  decrease  in  the  maximum  growth  rate  with  increased  adverse  shear  as  well  as  the 
increase  in  the  nondimensional  wavenumber,  ky0,  of  the  most  unstable  mode.  These  results 
agree  with  the  /3-plane  analyses  performed  numerically  by  Kuo  (1973)  and  analytically  by 
Schubert  et  a/.  (1992).  For  the  case  of  our  present  model,  A  will  remain  less  than  ~0.2. 

In  the  next  subsection,  we  integrate  the  shallow  water  equations  starting  with  various 
initial  vorticity  patterns  which  are  meant  to  represent  both  the  PV  strips  described  above 
and  the  ITCZ. 

4.3.2  Numerical  Experiments 

The  breakdown  of  PV  strips  for  nondivergent  flows  have  been  extensively  studied 
using  contour  dynamical  methods  (Dritschel,  1989;  Pratt  and  Pedlosky,  1991)  in  which 
the  positions  of  boundaries  between  between  regions  of  piecewise  constant  PV  are  pre¬ 
dicted.  Although  contour  dynamics  provides  an  elegant  method  for  studying  instabilities, 
it  requires  discontinuous  vorticity  patterns.  In  addition,  the  contour  dynamical  studies 
mentioned  above  were  both  balanced  studies.  In  this  subsection,  we  relax  the  assump¬ 
tion  of  balance  and  integrate  the  divergent  barotropic  (shallow  water)  equations  for  five 
different  experiments. 
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Figure  4.3:  The  change  in  normalized  growth  rates  and  nondimensional  wavenumber  with 
increased  adverse  shear.  The  solid  line  corresponds  to  normalized  growth  rates  which  are 
read  on  the  left  ordinate  and  the  dashed  line  corresponds  to  the  most  unstable  (nondi¬ 
mensional)  wavenumbers  which  are  read  on  the  right  ordinate. 
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Unlike  contour  dynamical  models,  our  model  uses  spectral  discretization.  This  means 
a  discontinuous  initial  vorticity  field  is  not  permissible  since  that  would  require  an  infinite 
number  of  spectral  coefficients  to  represent  the  jump.  We  therefore  represent  a  strip  of  PV 
in  two  different  ways.  The  first  method  defines  the  PV  strip  to  be  Gaussian  in  y  (centered 
about  the  line  y  =  yc)  and  constant  in  x,  i.e. 

C(y)  =  C  exp  (“ ,  (4-54) 

where  is  the  Gaussian  half-width  which  defines  the  e-folding  distance  and  £  is  the 
maximum  value  of  vorticity.  The  second  vorticity  pattern  to  be  defined  is  also  constant  in 
x  but  rather  than  being  Gaussian  in  y,  we  use  the  profile  function  (3.19)  to  describe  the 
strip’s  variation  in  y.  Once  the  vorticity  pattern  has  been  defined,  the  balanced  wind  and 
mass  fields  are  determined  by  the  nonlinear  balance  equation  described  in  Chapter  2.  Even 
though  the  initial  vorticity  pattern  is  unstable,  instabilities  will  not  grow  unless  the  fluid 
is  perturbed.  To  perturb  the  fluid  without  exciting  large  gravity  waves,  we  add  white 
noise  to  half  the  slow  mode  spectral  coefficients  and  leave  the  fast  modes  undisturbed. 
The  remaining  half  of  the  slow  mode  spectral  coefficients  are  computed  as  the  complex 
conjugates  of  the  disturbed  half.  This  ensures  the  physical  fields  will  remain  real. 

As  mentioned  previously,  the  mean  vorticity  field  must  first  be  removed  before  it  can 
used.  This  can  easily  be  done  by  transforming  the  vorticity  field  to  spectral  space,  setting 
the  spectral  coefficient  corresponding  to  k  =  l  =  0  to  zero,  and  transforming  the  field  back 
to  physical  space. 

Experiment  one  is  intended  to  simulate  a  continuous  ITCZ  pattern  which  is  centered 
at  10°N  in  the  eastern  Pacific  Ocean.  We  start  from  an  initial  state  which  is  defined  by 
a  PV  strip  extending  across  the  entire  width  of  the  domain.  The  strip  is  Gaussian  in  y 
with  a  half-width  of  235  km  and  a  maximum  vorticity  of  6.25  x  10“5  s-1(2.5/),  which 
gives  an  adverse  shear  of  ~0.065.  Figure  4.4  shows  the  height,  wind,  normalized  PV,  and 
absolute  vorticity  fields  at  the  initial  time.  The  PV  is  normalized  by  the  mean  depth  H 
which  is  taken  to  be  300  m.  This  value  of  H  provides  a  pure  gravity  wave  speed  of  ~54 
ms-1  and  a  Rossby  radius  of  ~1360  km  if  /  is  calculated  at  10°N.  The  absolute  vorticity 


Figure  4.4:  Dynamic  fields  at  initial  time  for  experiment  one:  height  (upper  left),  winds 
(upper  right),  normalized  PV  (lower  left),  and  absolute  vorticity  (lower  right).  The  height 
field  has  a  contour  interval  of  5  m.  The  maximum  wind  vector,  defined  as  the  length 
between  two  successive  tick  marks r  is  20  ms-1.  The  normalized  PV  and  absolute  vorticity 
have  contour  intervals  of  25xl0-8  and  are  in  units  of  10-8  s-1.  Tick  marks  are  every  400 
km. 
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(lower  right)  and  normalized  PV  fields  (lower  left)  have  contour  intervals  of  2.5  xlO-5  s-1 
or  1/  starting  at  2/.  These  contours  are  scaled  by  106.  The  height  field  (upper  right)  has 
contour  intervals  of  5  m  and  the  maximum  wind  vector  length  (as  measured  between  two 
vectors  without  allowing  overlap)  is  20  ms-1.  The  grid  size  is  6400  km  x 6400  km  with 
128x128  points  on  the  transform  grid,  thus  42  waves  are  kept  in  both  x  and  y.  The  initial 
wind  varies  from  —11  ms-1  to  11  ms-1  giving  a  total  wind  shear  across  the  ITCZ  oi‘  22 
ms-1.  Figures  4.5-4.8  show  the  evolution  of  the  vorticity  and  PV  fields  at  intervals  of  two 
days  starting  at  the  eighth  model  day.  All  fields  are  shown  for  day  fourteen  in  Fig.  4.9. 
Prior  to  the  eighth  day,  instabilities  are  not  easily  recognizable.  After  the  eighth  day,  the 
disturbances  grow  quite  rapidly.  Especially  interesting  to  note  is  the  development  of  two 
disturbances  placed  relatively  close  to  one  another.  As  time  proceeds,  these  disturbances 
undergo  a  merger  process  which  appears  to  be  complete  by  the  fourteenth  model  day. 
The  merger  of  geostrophic  vortices  in  laboratory  experiments  has  been  well  documented 
by  Griffiths  and  Hopfinger  (1987)  who  noted  that  vortices  will  merge  if  they  are  less  than 
3.3±0.2  radii  apart.  This  criterion  appears  to  be  satisfied  in  the  model  simulation.  By 
day  eleven  (not  shown),  the  two  closest  vortices  have  begun  the  merger  process.  By  day 
fourteen,  the  process  is  complete  and  only  two  of  the  three  original  disturbances  remain. 
The  origin  of  the  outer  bands  is  especially  noticeable  on  the  twelfth  day.  On  this  day, 
the  pooling  process  has  created  a  rather  thin  filament  of  PV  which  connects  the  circular 
vortex  to  the  merging  vortex.  This  filament  is  precisely  what  we  mean  by  an  outer  band. 
This  band  is  still  clearly  visible  in  the  PV  field  by  the  fourteenth  day. 

Before  continuing,  we  note  that  the  time  scale  of  the  model  ITCZ  breakdown  appears 
slower  than  that  observed  in  nature.  The  model  took  14  days  for  the  ITCZ  to  breakdown 
whereas  in  nature,  the  ITCZ  generally  breaks  down  in  three  to  seven  days.  It  must  be 
made  clear,  however,  that  in  nature  the  ITCZ  is  not  perturbed  with  such  small  amplitude 
forcings.  If  one  considers  day  eight  as  representative  of  a  finitely  perturbed  ITCZ  (i.e. 
letting  day  eight  represent  day  one  in  nature)  then  the  simulated  ITCZ  breakdown  is  very 
close  to  that  in  nature.  Faster  breakdowns  can  also  be  obtained  simply  by  increasing  the 
amplitude  of  the  perturbations. 


Figure  4.5:  Vorticity  (top)  anti  FV  (bottom)  at  day  eight  for  experiment  one.  Two  domain 
periods  are  shown.  Units  are  10-6  s-1  and  contour  intervals  are  25xl0~6. 
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Figure  4.6:  Same  as  for  Fig.  4.5  but  for  day  ten. 
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Figure  4.8:  Same  a8  for  Fig.  4.5  but  for  day  fourteen. 


Figure  4.9:  Same  as 
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Fig.  4.4  out  for  day  fourteen. 
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Figure  4.10:  Initial  relative  vorticity  profiles  for  experiments  one  through  three.  The 
smooth  line  is  for  Gaussian  profile  and  dashed  is  for  the  Melander  profile.  The  relative 
vorticity  has  units  of  s-1  and  is  scaled  by  105.  The  distance  is  measured  in  hundreds  of 
kilometers. 

For  the  second  experiment,  we  consider  the  effect  of  changing  the  shape  of  the  initial 
PV  profile.  We  choose  an  initial  vorticity  field  which  makes  use  of  the  profile  functions 
(3.19)-(3.20)  (hereafter  referred  to  as  a  Melander  profile)  and  has  the  same  average  vor¬ 
ticity  and  initial  winds  as  the  previous  Gaussian  experiment.  In  order  to  use  (3.19)  in 
Cartesian  coordinates,  we  simply  replace  r  by  |y  —  yc|.  The  profile  obtained  by  letting 
Ri  =  100  km  and  Rq  =  400  km  is  shown  in  Fig.  4.10  along  with  the  Gaussian  curve  used 
in  experiment  one.  By  examining  Fig.  4.10,  it  is  obvious  that  we  are  more  closely  approx¬ 
imating  the  uniform  vorticity  strips  which  were  considered  in  the  normal  mode  stability 
analysis. 

As  with  experiment  one,  we  again  perturb  the  above  vorticity  distribution  by  adding 
white  noise  to  the  balanced  modes.  It  is  important  to  note  that  the  same  sequence  of 
random  numbers  is  used  for  both  experiments  one  and  two.  The  effect  of  changing  this 
will  be  examined  in  the  final  experiment.  Results  from  experiment  two  for  days  ten  and 
fourteen  are  shown  in  Figs.  4.11-4.13.  Both  experiments  show  the  same  trend  of  ITCZ 
breakdown  and  the  formation  of  outer  bands.  The  primary  difference  between  experiments 
one  and  two  is  the  observed  appearance  of  a  third  vortex.  Although  the  third  vortex  did 


Figure  4.11:  Vorticity  (top)  and  PV  (bottom)  at  10  days  for  experiment  two.  Units  and 
contour  intervals  are  the  same  as  for  experiment  one. 
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Figure  4.12:  Same  as  Fig.  4.11  but  for  day  fourteen. 


Figure  4.13:  Dynamic  fields  for  experiment  two  at  fourteen  days.  Panels,  units,  and 
contour  intervals  are  the  same  as  in  experiment  one. 
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not  appear  in  experiment  two,  it  is  interesting  to  note  that  the  areal  extent  of  the  pooled 
vorticity  is  greater  near  the  location  where  the  two  merging  vortices  formed  in  experiment 
one.  It  is  also  clear  that  outer  bands  form  regardless  of  the  initial  shape  of  the  vorticity 
field. 

For  the  third  experiment,  we  examine  the  effect  of  increasing  the  width  of  the  PV 
strip  on  both  the  breakdown  of  the  ITCZ  and  the  formation  of  outer  bands.  The  initial 
vorticity  field  for  experiment  three  is  given  as  a  Gaussian  distribution  in  y  with  a  half¬ 
width  of  300  km  and  a  maximum  vorticity,  £,  of  2/.  This  value  of  £  again  produces  a 
maximum  zonal  wind  of  ~11  ms-1.  Again,  as  in  experiments  one  and  two,  we  perturb 
the  balanced  modes  with  the  same  sequence  of  random  numbers.  The  model  results 
are  shown  in  Figs.  4.14-4.16.  As  expected  from  linear  theory,  the  weaker  C  produces 
slower  growth  rates  with  experiment  three  taking  sixteen  model  days  for  the  vortices  to 
fully  detach  as  opposed  to  fourteen  days  for  the  earlier  experiments.  Again  from  linear 
theory,  we  should  also  expect  wider  vorticity  distributions  to  produce  greater  disturbance 
wavelengths.  Comparison  with  experiments  one  and  two,  however,  shows  no  noticeable 
increase  the  the  wavelengths  of  the  disturbances.  This  may  partially  be  explained  by  the 
model’s  boundary  conditions  which  force  periodicity.  Additionally,  the  wavelength  of  the 
most  unstable  mode  does  not  vary  by  a  significant  amount  between  experiments  one  and 
three.  This  can  be  seen  in  Fig.  4.17  which  gives  the  nondimensional  growth  rate  as  a 
function  of  the  zonal  wavenumber  s,  where  s  is  defined  as  the  domain  length  in  the  i 
direction  divided  by  the  wavelength.  Thus,  s  gives  the  maximum  number  of  waves  for 
a  given  wavelength  that  can  fit  in  the  models  x  domain.  Figure  4.17  shows  that  the 
most  unstable  wavenumber  varies  little  between  experiments  one  and  three.  It  therefore 
seems  reasonable  that  because  the  boundary  condition  forces  periodicity,  both  the  y/,= 235 
km  and  the  y/,=300  km  Gaussian  strips  could  easily  produce  wavenumber  two  patterns. 
However,  as  the  strips  are  made  thinner,  they  are  more  inclined  to  produce  wavenumber 
three  type  patterns.  In  making  these  statements,  we  must  remember  that  comparison  of 
nonlinear  results  with  linear  normal  mode  theory  is  highly  subjective.  The  half-width  of 
a  Gaussian  profile  may  not  adequately  approximate  the  half-width  of  a  uniform  vorticity 
profile. 
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Figure  4.14:  Vorticity  (top)  and  PV  (bottom)  at  twelve  days  for  experiment  three.  Units 
and  contour  intervals  are  the  same  as  for  experiments  one  and  two. 
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Figure  4.15:  Same  as  Fig.  4.14  but  for  day  sixteen. 
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Figure  4.16:  Dynamic  fields  for  experiment  three  at  sixteen  days.  Panels,  units,  and 
contour  intervals  are  the  same  as  in  experiment  one  and  two. 
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Figure  4.17:  Nondimensional  growth  rate  versus  nondimensional  zonal  wavenumber  s  (see 
text  for  detail).  Dotted  line  corresponds  to  yo  —  300  km,  dashed  line  corresponds  to 
yo  =  235  km,  and  solid  line  corresponds  to  yo  =  200. 

In  experiment  four  we  test  the  sensitivity  of  experiment  one’s  results  to  the  initial 
condition.  We  do  this  by  simply  changing  the  seed  of  the  random  number  generator,  thus 
producing  a  new  sequence  of  random  numbers,  and  integrating  the  model  as  before.  The 
results  of  this  experiment  for  days  ten  and  fourteen  are  presented  in  Fig.  4.18  and  Fig. 
4.19.  These  results  demonstrate  the  tremendous  sensitivity  of  the  model  to  the  initial 
condition.  The  linear  analysis  can  again  help  us  to  understand  these  results.  As  we  noted 
in  experiment  one,  Fig.  4.17  shows  that  growth  rates  between  two  wavenumbers  vary  little, 
especially  for  thinner  strips.  Thus,  if  the  perturbation  excites  wavenumber  two  slightly 
more  than  wave  number  three,  a  wave  number  two  pattern  would  be  expected  and  vice 
versa.  Figure  4.17  also  suggests  that  using  the  Gaussian  half-width  as  an  approximation 
to  a  uniform  vorticity  fields  half-width  may  not  be  accurate.  Results  from  the  nonlinear 
integration  suggest  that  a  width  considerably  shorter  than  the  Gaussian  half- width  may 
provide  a  better  approximation. 

In  the  previous  experiments  we  assumed  the  ITCZ  was  a  continuous  disturbance 
across  the  entire  domain  of  the  model.  We  now  wish  to  consider  the  effects  of  a  localized 
ITCZ.  Therefore,  for  the  fifth  experiment,  we  consider  an  initial  vorticity  field  which  is 
Gaussian  in  y  (with  y h  =  235  km)  but  confined  in  x.  To  accomplish  this,  we  again  make 
use  of  the  Melander  profile  equations.  We  now  specify  the  profile  function  to  be  function 
of  |x  —  xc|  and  let  Rq  =  1600  km  and  Ri  =  2000  km.  This  produces  a  disturbance  which 


116 


Figure  4.18:  PV  fields  for  experiment  four  at  10  days.  The  three  PV  strips  were  each 
generated  using  the  same  initial  vorticity  field,  but  each  were  perturbed  using  different 
sequences  of  random  numbers  with  identical  ranges  in  magnitude.  The  middle  PV  strip 
corresponds  to  experiment  one.  The  domain  is  plotted  twice.  Units  and  contour  intervals 
are  the  same  as  for  experiment  one. 
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Figure  4.19:  Same  as  Fig.  4.18  but  at  fourteen  days. 
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is  approximately  one  half  of  the  domain  length.  The  maximum  value  of  vorticity  is  again 
chosen  to  be  2.5/.  The  model  is  run  for  a  period  of  four  days  by  which  time  the  strip 
has  broken  down  into  two  separate  vortices.  The  results  of  the  four  day  calculation  are 
shown  in  Figs.  4.20-4.23.  Again,  as  in  experiments  one  through  four,  the  ITCZ  breaks 
down  and  thin  bands  of  PV  form  between  the  incipient  vortices.  By  day  four,  all  fields 
indicate  the  existence  of  two  separate  vortices.  This  is  roughly  four  times  faster  than 
the  experiments  with  a  continuous  disturbance.  This  is  not  surprising,  however,  since 
wave  breaking  at  large  x  is  aiding  the  ITCZ  breakdown  process  by  providing  rather  large 
perturbations  in  the  PV  field.  Another  noticeable  difference  between  experiments  five  and 
the  others  is  the  motion  of  the  disturbances.  Clearly,  the  circulation  of  each  incipient 
vortex  is  influencing  the  other’s  motion;  thus,  the  two  vortices  simply  rotate  around  one 
another  (the  Fujiwhara  effect).  This  same  feature  appears  to  be  evident  in  the  satellite 
images  shown  in  Figs.  4. 1-4.2  (although,  admittedly,  spherical  effects  are  also  largely 
responsible  for  the  vortex  positions). 

To  summarize  our  results,  PV  dynamics  and  barotropic  instability  theory  provide 
useful  tools  for  examining  both  the  breakdown  of  the  ITCZ  and  the  formation  of  outer 
bands.  This  in  no  way  implies  other  concepts,  such  as  moist  dynamical  processes  and 
baroclinic  effects,  are  not  important.  Clearly,  they  must  be.  We  have  simply  isolated  a 
fundamental  part  of  the  dynamics  and  have  demonstrated  its  usefulness  in  understanding 
the  complex  processes  of  ITCZ  breakdown  and  outer  band  formation. 

In  the  next  chapter,  we  consider  the  stability  of  spiral  bands  after  their  formation. 


Figure  4.20:  Dynamic  fields  at  initial  time  for  experiment  five.  Panels,  units,  and  contours 
are  identical  to  experiments  one  through  four. 
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Figure  4.22:  Same  as  figure  4.20  but  for  day  three 


Figure  4.23:  Same  as  figure  4.20  but  for  day  four. 


Chapter  5 


ON  THE  STABILITY  OF  SPIRAL  BANDS 

In  Chapter  3  we  described  the  evolution  of  spiral  bands  in  tropical  cyclones  and 
noted  the  ejection  of  PV  tongues  into  areas  of  relatively  low  PV.  In  Chapter  4,  we  noted 
that  reversed  PV  gradients  satisfied  the  necessary  condition  for  barotropic  instability  and 
helped  explain  the  breakdown  of  the  ITCZ  into  pools  of  PV.  This  raises  the  the  following 
question.  Why  are  spiral  bands  are  never  observed  to  be  unstable  even  though  they, 
themselves,  are  most  probably  areas  of  reversed  PV  gradients.  In  this  chapter  we  attempt 
to  shed  some  light  on  the  this  problem  by  examining  the  stability  of  PV  rings  surrounding 
a  main  vortex.  We  attempt  to  show  that  the  main  vortex  of  the  tropical  cyclone  provides 
enough  adverse  shear  to  resist  the  “rolling  up”  of  spiral  bands.  We  also  examine  the 
stability  of  annular  regions  of  high  PV  in  the  absence  of  any  main  vortex.  This  will 
help  shed  some  light  on  such  observed  features  as  hurricane  meso-vortices  and  polygonal 
eyewalls. 

We  begin  by  considering  the  normal  mode  stability  analysis  of  multiple  concentric 
annular  regions  of  PV.  In  section  5.2  we  consider  the  specific  case  of  three  regions  of 
different  PV.  This  case  involves  a  center  region  of  irrotational  air,  a  ring  of  high  PV  air, 
and  a  surrounding  region  of  irrotational  air  as  r  goes  to  infinity.  We  use  this  case  to 
examine  the  stability  of  bands  without  the  presence  of  a  main  vortex.  This  case  also 
has  dynamic  links  to  eyewall  instabilities  and  polygonal  eyewalls.  In  the  final  section,  we 
consider  a  four  region  model.  The  fourth  region  is  the  addition  of  a  high  PV  region  in  the 
center  of  the  vortex.  We  use  this  to  test  the  effects  of  intense  adverse  shear  on  the  growth 
rates  of  instabilities. 
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Figure  5.1:  The  piecewise  continuous  mean  vorticity  distribution  as  a  function  of  r.  As 
one  proceeds  radially  inwards,  the  vorticity  jumps  in  increments  of  . 

5.1  Normal  Mode  Stability  Analysis 

In  this  section  we  derive  a  method  for  determining  the  normal  mode  stability  of  annu¬ 
lar  PV  strips  to  small  amplitude  disturbances.  We  consider  the  basic  state  axisymmetric 
vorticity  distribution  shown  in  Fig.  5.1.  At  r  =  rj  (j  =  0, 1, 2, . . . ,  J)  the  vorticity  jumps 
by  an  amount  (j  as  one  proceeds  radially  inward,  where  £j  can  be  either  positive  or  neg¬ 
ative.  The  basic  state  streamfunction  corresponding  to  this  vorticity  distribution  is  given 
by 

j  J 

$(r)  =  £  ln(r/ry)  +  £  (r2  -  rfj  for  ry  <  r  <  rj+l,  (5.1) 

j'=l  j'-j+i 

which  can  be  confirmed  by  two  differentiations,  the  first  of  which  yields 

j  J 

™(r)  =  H  k>'r2  for  rj  <r^  Ti+u  (5.2) 

j'=i  j'=j+i 

and  the  second  of  which  yields 

J 

C(r)  =  H  Zj’  for  Tj<r<  rJ+i. 

i'=i+ 1 


(5.3) 
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Since  (5.3)  corresponds  to  the  stairstep  vorticity  pattern  shown  in  Fig.  5.1,  and  since  rj>(r) 
and  t7(r)  are  continuous,  (5.1)  is  the  solution  of  the  axisymmetric  barotropic  invertibility 
principle  d/rdr(rdrp/dr)  =  £. 

Now  suppose  that  each  interface  is  perturbed  by  an  amount  r)j(<p,t)  in  the  sinusoidal 
fashion  rjj(<p,t)  =  The  solution  for  the  perturbation  streamfunction  is 

i  J 


13  Ai'  (ri'/r)m  +  S  Ai'irlri,y 

b'=i  i'=j+i 


eHm<p-i/t)  rj  <r  <  Tj+i.  (5.4) 


It  can  easily  be  confirmed  that  t /)'  is  continuous.  Continuity  of  v  at  r_,  +  r]j  yields  mAj  = 
-fojilj.  Using  this  in  (5.4),  and  the  resulting  expression  in 


dr\j  _  dru  im  , 

+u‘at  =  -—' *{ri)' 


dt 


we  obtain 


where 


(u  -  mu>j)fij  +  23  =  0, 

j’=i 


(5.6) 


(5.7) 


j(m)  =  (  (rr/ri)  r<3 

**  1  (rj/rj')ml  f>3 , 

and  Qj  s  v/rj.  This  form  of  the  eigenvalue/eigenvector  problem  is  identical  to  that 

presented  by  Dritschel  (1989). 

If  all  the  £ji  vanish  except  the  one  for  j'  =  j,  we  obtain 

m  1  2m  2  1  \  mj 

which  is  the  special  case  found  in  our  discussion  of  PV  waves  in  Chapter  3. 

5.2  The  Stability  of  Annular  Potential  Vorticity  Regions 

We  now  consider  the  specific  case  of  a  single  annular  region  of  PV  (i.e.  J  =  2),  there 
are  two  interfaces  for  this  case  as  shown  in  Fig.  5.2.  The  eigenvalue  problem  for  this 
simplified  case  reduces  to 

™ i  -  36  -^2(n/r2)m-1  \  /  m  \  _  (  m 

i(nA2)m+1  mu>2  —  3^2  /  \  V2  )  \  m 


( 


(5.8) 


Figure  5.2:  The  assumed  vorticity  profile  for  J  —  2. 


If  we  note  that  a>i  =  0,  u>2  =  j[(ri/r2)2  —  1],  and  =  —£2  for  this  special  case,  it  can  easily 
be  verified  that  u  will  be  pure  real  for  m  =  0, 1,  or  2.  This  implies  the  vorticity  field  will 
remain  stable  to  these  disturbance  patterns.  The  remaining  wavenumbers  can,  however, 
produce  frequencies  with  imaginary  components.  Figure  5.3  shows  the  normalized  growth 
rates  1 /,-/&  as  a  function  of  the  nondimensional  width  of  the  annular  region,  n/r2.  Clearly, 
thinner  annular  regions  should  produce  the  highest  growth  rates  but  at  much  higher 
wavenumbers.  The  nonlinear  normal  mode  analysis  also  provides  information  on  the 
amplitude  and  phase  of  the  most  unstable  modes.  Fig.  5.4  shows  the  most  unstable 
amplitude  and  phase  of  a  wave  number  four  pattern  with  ri  =  200  km  and  r2  =  250  km. 

To  test  our  results  numerically  for  the  divergent  barotropic  case,  we  initialize  our 
model  with  an  annular  region  which  is  Gaussian  in  r.  Unlike  in  Chapter  4,  however,  we 
cannot  simply  perturb  the  model  with  white  noise.  The  reason  for  this  is  that  experience 
shows  that  the  influence  of  the  doubly  periodic  boundary  condition  will  always  force  the 
growth  of  a  wavenumber  four  pattern.  To  overcome  this  problem,  we  instead  perturb  the 
model  with  a  weak  forcing  of  the  wavenumber  desired.  The  wavenumber  of  the  forcing 
is  chosen  to  be  the  most  unstable  wavenumber  based  on  linear  theory  while  using  the 
the  Gaussian  half-width  to  represent  the  width  of  the  ring.  Although,  we  cannot  say  for 
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Figure  5.3:  The  normalized  growth  rates  as  a  function  of  the  width  of  the  annular  region 
for  wavenumber  3,  4,  5,  and  6  disturbances. 

certain  that  the  chosen  wavenumber  is  necessarily  the  most  unstable  for  the  nonlinear 
divergent  case,  it  does  allow  us  to  view  the  evolution  of  various  unstable  patterns. 

To  demonstrate  the  growth  of  instabilities  within  an  annular  region,  we  integrate  the 
model  starting  from  an  initial  state  has  zero  vorticity  everywhere  except  within  an  annular 
region.  Within  the  annular  region,  the  vorticity  is  Gaussian  in  r  with  a  half-width  of  35 
km  and  a  maximum  value  of  2.5xl0-4  s-1.  The  center  of  Gaussian  profile  is  located  at 
a  radius  of  225  km.  If  the  Gaussian  half-width  is  taken  to  be  representative  of  the  half¬ 
width  a  uniform  annular  region,  then  the  most  unstable  wavenumber  for  this  case  should 
be  four,  i.e.  t\/t2  —  0.73.  For  this  simulation  and  all  others  in  the  chapter,  we  used  1282 
points  on  the  transform  grid  and  kept  42  waves.  The  initial  fields  for  this  experiment  are 
shown  in  Fig.  5.5  to  demonstrate  the  amplitude  of  the  initial  forcing.  The  results  from 
this  experiment  are  shown  in  Figs.  5.6-5.8  for  24,  36  and  48  hrs  respectively.  By  24 
hours  the  growing  instabilities  are  clearly  visible  in  both  the  PV  and  vorticity  fields.  The 
patterns  in  both  these  fields  bear  remarkable  similarity  to  the  pattern  predicted  by  linear 
theory  (Fig.  5.4).  By  36  hours,  four  distinct  circulation  patterns  have  clearly  developed 
as  evidenced  by  all  fields.  These  four  vortices  are  completely  detached  by  48  hours. 

Similar  results  were  also  found  for  wavenumber  three  and  five  perturbations.  These 
results  are  shown  in  Figs.  5.8-5.11.  For  the  wavenumber  three  perturbation,  we  used  a 


Figure  5.4:  The  amplitude  and  phase  relationship  for  a  wavenumber  four  disturbance  in 
which  ri  =  200  km  and  ri  —  250  km.  The  amplitudes  of  the  disturbances  are  known  only 
to  within  an  arbitrary  constant. 

Gaussian  half-width  of  50  km  centered  at  a  radius  of  225  km.  Again  if  assume  the  Gaussian 
half-width  is  representative  of  the  half-width  of  a  uniform  annular  region,  then  we  obtain 
a  ratio  of  t\/t2  =  0.63,  which  is  most  unstable  for  wavenumber  three  disturbances.  The 
wavenumber  five  perturbation  used  a  Gaussian  half-width  of  20  km  centered  at  225  km 
which  gives  a  ratio  of  ri/r2  =  0.83.  The  effect  of  increasing  the  width  of  the  disturbance 
follows  linear  theory  well.  As  the  width  of  the  disturbance  is  increased  the  growth  rates 
of  the  instabilities  decrease.  This  is  easily  confirmed  by  comparing  Fig.  5.9  with  Fig.  5.11, 
and  Fig.  5.10  with  Fig.  5.12.  Clearly,  the  thinner  strip  evolved  faster. 

Although  we  have  been  examining  the  stability  of  these  annular  regions  with  regards 
to  spiral  bands,  the  results  may  also  have  applications  to  the  eyewall  of  a  tropical  cyclone. 
The  vorticity  field  for  the  eyewall  of  a  tropical  cyclone  is  similar  in  shape  to  the  above 
annular  regions  of  vorticity  in  that  the  peak  vorticity  is  not  observed  in  the  center  of  the 
vortex  (as  we  have  assumed  in  the  previous  chapters)  but  at  some  small  distance  away 
from  the  center.  It  is  therefore  possible  for  instabilities  of  the  type  discussed  above  to  be 
found  near  the  eyewall.  There  are  two  types  of  observations  which  support  such  a  theory. 
The  first  is  the  observed  phenomena  of  meso-vortices  in  some  tropical  cyclone  eyewalls 
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Figure  5.5:  The  initial  wind  and  mass  fields  for  the  wavenumber  four  perturbed  case.  The 
height  field  (upper  left)  has  contour  intervals  of  2  m.  The  wind  field  (upper  right)  has  a 
maximum  wind  vector  of  12ms~l.  The  normalized  PV  (lower  left)  and  absolute  vorticity 
(lower  light)  are  in  units  of  10~5  s-1  and  have  contour  intervals  of  l.OxlO-4. 
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Figure  5.7:  The  same  as  Fig.  5.5  but  for  36  hours. 
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Figure  5.9:  The  wind  and  mass  fields  for  the  wavenumber  three  perturbed  case  at  24  hrs. 
The  maximum  wind  vector  is  15  ms-1.  Panels,  units  and  contour  intervals  are  the  same 
as  for  Fig.  5.5. 
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Figure  5.10:  The  same  as 
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'ig.  5.9  but  for  48  hours. 
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Figure  5.11:  The  wind  and  mass  fields  for  the  wavenumber  five  perturbed  case  at  24  hrs. 
The  maximum  wind  vector  is  8  ms"1.  Panels,  units,  and  contour  intervals  are  the  same 
as  for  t  lg.  5.5. 
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Figure  5.12:  The  same  as  Fig.  5.11  but  for  48  hours. 
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(Black  and  Marks,  1991).  The  second  is  the  polygonal  eyewall  shapes  observed  in  some 
tropical  cyclones  (Lewis  and  Hawkins,  1982;  Muramatsu,  1986). 

Black  and  Marks  (1991)  observed  a  meso- vortex  in  hurricane  Hugo  (1989)  during  what 
was  intended  to  be  a  routine  penetration  at  450  m.  The  vortex  had  wind  and  pressure 
perturbations  of  ~30ms_1  and  12  mb,  respectively.  The  perturbation  winds  were  roughly 
a  factor  of  three  smaller  than  than  those  in  the  eyewall.  The  perturbation  vorticity, 
however,  was  observed  to  be  about  one  order  of  magnitude  greater  and  nearly  equivalent 
to  that  of  a  tornado.  The  above  model  results  also  indicate  similar  trends.  That  is, 
the  vorticity  is  highest  in  the  model  vortices  even  though  the  winds  are  not  significantly 
greater  and  the  geopotential  is  also  clearly  less  in  the  vortices.  Although  the  dynamics  of 
a  true  meso- vortex  are  considerably  more  complicated  than  our  simple  model  can  handle, 
the  similarities  are  curious. 

Support  for  our  instability  theory  also  comes  from  the  observed  phenomena  of  polyg¬ 
onal  eyewalls.  As  evidenced  in  the  above  simulations,  an  annular  region  of  PV  can  evolve 
into  polygonal  patterns.  Lewis  and  Hawkins  (1982)  showed  there  to  be  distinct  polygonal 
patterns  in  the  eyewalls  of  several  hurricanes.  Figure  5.13  shows  an  example  of  two  such 
polygonal  eyewall  features  of  Hurricane  Betsy  (1965)  as  seen  by  WSR-57  radar.  Polygonal 
patterns  were  also  found  in  Hurricane  Anita  (1977),  Hurricane  Debbie  (1969)  and  Hurri¬ 
cane  Frederic  (1979).  To  prove  these  polygonal  features  were  not  simply  artifacts  of  the 
radar,  Lewis  and  Hawkins  (1982)  were  able  to  obtain  simultaneous  data  from  both  the 
National  Weather  Service’s  linear  receiver  radar  and  the  RFC  airborne  digital  log  receiver 
radar  system  for  Hurricane  Anita.  Both  systems  showed  similar  pentagon  eye  patterns. 
In  addition,  Lewis  and  Hawkins  showed  that  the  polygonal  patterns  also  appeared  in 
the  time-integrated  radar  rainfall  estimates  of  Hurricane  Frederic,  again  suggesting  that 
these  polygonal  eyewall  patterns  are  real  physical  features  and  not  simply  artifacts  of  the 
observing  system. 

The  dynamical  theory  put  forth  by  Lewis  and  Hawkins  was  that  these  patterns  were 
due  to  the  modulation  of  convection  by  horizontally  propagating  internal  gravity  waves. 
They  proposed  that  interference  patterns  due  to  the  superposition  of  differing  wavenum¬ 
bers  and  periods  would  tend  to  produce  polygonal  eyes.  Their  model  results,  however, 


Figure  5.13:  Two  examples  of  polygonal  eyewall  features  of  Hurricane  Betsy  as  seen  by 
WSR-57  radar.  (From  Lewis  and  Hawkins,  1982.) 

suggest  that  if  this  were  the  case,  polygonal  features  should  also  be  equally  prominent  in 
spiral  bands  as  well.  Their  observational  studies,  however,  show  only  one  such  case  where 
this  was  observed  (Caroline,  1975). 

In  a  similar  study,  Muramatsu  (1986)  performed  an  extensive  analysis  of  the  polygonal 
eyewall  patterns  in  typhoon  8019  WYNNE  (1960).  Muramatsu  also  performed  an  in- 
depth  comparison  of  several  hurricanes  and  typhoons  which  exhibited  polygonal  eyewall 
features,  the  results  of  which  are  shown  in  Table  5.1.  Muramatsu  determined  the  most 
common  characteristics  displayed  by  these  tropical  cyclones  were:  1)  concentric  eyes  in 
the  developed  tropical  cyclone  stage,  2)  central  pressures  ranging  from  920-950  mb.  and 
3)  square  to  hexagonal  shapes  occurred  in  highest  frequency  (pentagon  was  observed 
most  frequently  and  persisted  the  longest)  .  These  observations  lend  support  to  the 
instability  theory  presented  above.  Concentric  eyewalls  suggest  definite  peaks  in  the  PV 
field,  thus  more  closely  resembling  annular  PV  regions.  Secondly,  the  more  intense  tropical 
cyclones,  as  determined  by  central  pressure,  would  be  expected  to  have  higher  PV  values 


139 


Table  5.1:  Summary  of  tropical  cyclone  characteristics  and  features  of  polygonal  eyes. 
(Prom  Muramatsu,  1986). 


Cyclone  names 

Feature 

Eye  diameter 

inner  outer 

Concentric 
double  eye 

CeatraJ 

pressure 

DONNA  I960 

S 

23  km 

90  km 

ye* 

940  mb 

BETSY  1965 

6 

40 

120 

ye* 

950 

4 

40 

100 

ye* 

950 

DEBBIE  1969 

5 

20 

36 

ye* 

950 

ANITA  1977 

4 

33 

72 

ye* 

926 

S 

33 

67 

ye* 

926 

FREDERIC  1979 

4 

30 

100 

ye* 

950 

ALLEN  1980 

6 

24 

130 

ye* 

920 

Typhoon  8019 

4-6 

30 

260 

yet 

920-935 

WYNNE 

Typhoon  6118 

4-5 

26-35 

150 

ye* 

920 

and  therefore  be  more  susceptible  to  instabilities.  Lastly,  the  smaller  number  polygons 
should  be  observed  most  frequently  because  of  the  ratio  of  the  inner  to  outer  distance  of  the 
eyewalls.  However,  why  pentagons  should  be  favored  is  not  yet  understood.  A  weakness 
of  our  theory  is  that  Muramatsu  did  not  observe  any  tendency  of  the  outer  eyewalls  to 
also  show  polygonal  features,  although  Lewis  and  Hawkins  (1982)  did  show  one  such  case 
(Hurricane  Debbie,  1969).  The  numerical  experiments  we  performed  indicated  that  the 
vorticity  field  contained  minimal  gravity-wave  activity.  This  suggests  that  the  theory  of 
interfering  gravity  waves  seems  less  plausible  than  our  present  theory.  In  addition,  our 
present  theory  does  not  require  the  spiral  bands,  themselves,  to  be  polygonal  shaped. 
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Figure  5.14:  A  schematic  representation  of  adverse  shear.  The  shear  of  the  vortex  created 
by  the  decrease  in  tangential  wind  speed  with  radius  clearly  opposes  the  circulation  induced 
by  the  PV  anomaly. 

5.3  The  Effect  of  Adverse  Shear  on  Stability 

In  the  previous  section,  we  considered  the  growth  of  instabilities  for  a  annular  band 
which  was  free  of  influence  from  any  other  circulations.  Chapter  4  indicated  that  under 
these  situations,  the  PV  field  is  unstable  and  thus  instabilities  should  freely  grow.  In¬ 
deed,  this  was  the  case.  In  nature,  however,  the  spiral  bands  are  influenced  greatly  by 
circulation  of  the  main  vortex.  Since  the  tangential  winds  of  the  main  vortex  decrease 
with  radius  (beyond  the  radius  of  maximum  winds),  they  introduce  a  shear  which  opposes 
the  cyclonic  circulation  induced  by  the  PV  anomaly  of  the  spiral  band.  This  situation  is 
shown  schematically  in  Fig.  5.14.  When  the  shear  environment  opposes  the  anomaly,  we 
refer  to  this  as  adverse  shear  (Dritschel,  1989).  We  will  demonstrate  in  this  section  that 
the  tendency  for  instabilities  to  grow  in  spiral  bands  is  severely  hampered  by  adverse  shear 
of  the  vortex.  Our  discussion  begins  again  with  a  look  at  normal  mode  stability  analysis 
and  is  followed  by  a  numerical  simulation. 

5.3.1  Normal  Mode  Stability  Analysis 

We  simulate  the  effect  of  a  main  vortex  on  a  spiral  band  by  considering  an  annular 
region  of  vorticity  (as  before)  with  another  region  of  higher  vorticity  located  at  the  center 
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Figure  5.15:  An  example  of  the  vorticity  profile  for  J  —  3. 


of  the  annulus  (see  Fig.  5.15).  This  is  done  by  considering  the  case  of  J  —  3  in  (5.7),  for 
which  the  matrix  eigenvalue  problem  can  be  written 

f  -J^ri/rj)—1  -^6(nA3)m_1\  { Vl\  /»7i> 

-ift(n/*,a)m+l  -  ^2  -|6(»’2/>*3)m“1  V2  =  v  in  .  (5.9) 

k ~^l(rl/r3)m+l  -J6(fa/»'»)m+1  rnw 3-^3  )  \m)  \m> 

We  now  consider  the  specific  case  of  rj  =  50  km,  t?  =  200  km,  =  250  km,  £2  —  —£3  = 
2.5  x  10~4  8~l,  and  £1  is  allowed  to  vary  from  0.0  to  15.0  x  10-4  s-1.  The  frequencies 
and  corresponding  e-folding  times  of  the  most  unstable  modes  as  a  function  £1  are  shown 
in  Fig.  5.16.  When  £1  =  0.0,  the  results  are  identical  to  case  of  J  =  2  for  the  same 
values  of  ri  and  ri.  As  the  vorticity  of  the  main  vortex  increases,  the  adverse  shear  also 
increases.  Clearly,  the  growth  rates  of  the  most  unstable  modes  decrease  as  the  adverse 
shear  generated  by  the  main  vortex  increases.  Also  of  importance  is  the  trend  for  the 
wavenumber  of  the  most  unstable  mode  to  increase  as  the  adverse  shear  increases.  The 
implications  of  this  result  is  that  the  growth  rates  of  instabilities  in  spiral  bands  should 
be  expected  to  decrease  dramatically  as  the  strength  of  the  main  vortex  increases  and  the 
wavenumber  of  the  disturbance  should  increase.  The  fact  that  the  PV  in  spiral  bands  is 
never  observed  to  “pool”  supports  this  result. 
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Figure  5.16:  The  growth  rates  (u)  of  the  most  unstable  modes  for  wavenumbers  m  =3-9 
as  a  function  of  fr.  Constant  values  include  ri  =  50  km,  r2=200  km,  r$  =  250  km,  and 
£2  =  -u  =  2.5  x  10"4. 

5.3.2  Numerical  Simulation  of  Adverse  Shear 

To  demonstrate  the  effects  of  adverse  shear  on  a  spiral  band,  we  have  integrated  the 
wavenumber  four  case  of  the  previous  subsection  but  with  the  addition  of  an  initial  central 
vorticity  region  which  again  takes  the  form  of  the  Melander  profile.  The  Melander  profile 
for  this  case  is  defined  by  Ri  =  50  km,  Rq  =  100  km,  k  =  2.56,  and  (  =  5.0  x  10-4. 
Thus,  the  vorticity  within  the  main  vortex  is  only  twice  that  of  the  band.  To  keep 
comparisons  with  the  results  of  the  previous  section  fair,  we  have  also  given  the  model 
band  a  wavenumber  four  perturbation  initially.  Figure  5.17  shows  the  initial  fields  from 
this  experiment.  The  results  at  64  and  120  hrs  are  shown  in  Figs.  5.18  and  5.19.  Prior 
to  64  hours,  instabilities  in  the  band  are  not  noticeable.  Only  at  64  hours  and  beyond 
are  they  clearly  visible.  When  no  main  vortex  was  present,  the  wavenumber  four  pattern 
had  completely  detached  by  48  hours.  By  120  horns,  individual  pools  of  PV  are  clearly 
discernible  and  the  band  appears  to  display  a  wavenumber  eight  PV  pattern.  The  trend 
predicted  by  linear  normal  mode  theory  is  therefore  correct.  Not  only  did  the  instabilities 
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Figure  5.17:  The  initial  wind  and  mass  fields  for  the  wavenumber  four  perturbed  case  with 
a  main  vortex.  The  height  field  (upper  left)  has  contour  intervals  of  4  m.  The  wind  field 
(upper  right)  has  a  maximum  wind  vector  of  17ms-1.  The  PV  (lower  left)  and  vorticity 
(lower  right)  are  in  units  of  10-5  s-1  and  have  contour  intervals  of  l.OxlO-4. 
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Figure  5.18:  The  same  as  Fig.  5.17  but  for  64  hours. 
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take  longer  to  amplify,  but  much  higher  wavenumbers  were  produced.  In  nature  we  should 
expect  these  results  to  be  even  more  dramatic  because  of  the  extremely  high  vorticity 
associated  with  the  main  vortex  of  a  tropical  cyclone,  which  have  been  measured  to  be  as 
high  as  4.8  x  10~3  s-1  (Hawkins  and  Imbembo,  1976).  Although  numerous  calculations  of 
vorticity  within  spiral  bands  are  not  available,  results  from  Powell  (1990)  and  Ryan  et  a l. 
(1992)  suggest  the  vorticity  in  the  bands  can  range  from  1.5  xlO-4  s-1  to  2.5  xlO-3  s-1. 
If  these  values  are  assumed  to  representative  of  typical  hurricane  bands,  the  lifetime  of 
the  bands  should  be  much  less  than  the  time  required  for  instabilities  to  affect  the  spiral 
bands.  Even  if  the  vorticity  is  in  the  main  vortex  is  only  twice  that  of  the  band,  as  in  the 
above  experiment,  it  is  doubtful  that  a  single  band  could  persist  in  nature  for  120  hours 
(5  days)  as  required  by  model  results. 


Chapter  6 


SUMMARY  AND  CONCLUSIONS 

We  have  presented  several  new  theories  concerning  the  formation  and  stability  of  both 
inner  and  outer  hurricane  spiral  bands  through  use  of  a  limited  area,  shallow  water,  normal 
mode  spectral  model.  Because  the  basis  functions  of  the  model  were  the  normal  modes  of 
the  system,  we  were  able  to  naturally  partition  the  geostrophic  (rotational)  modes  from  the 
gravity-inertia  modes.  The  bands  were  found  to  project  almost  entirely  onto  the  rotational 
modes.  These  results  indicate  that  the  gravity  wave  theories  conventionally  accepted  for 
hurricane  bands  (e.g.  Tepper,  1958,  Abdullah,  1966;  Kurihara,  1976;  Willoughby,  1978) 
may  not  be  entirely  appropriate.  Evidence  was  supplied  to  suggest  that  spiral  bands 
should  more  appropriately  be  classified  as  slow  manifold  phenomena  (Leith,  1980;  Daley, 
1991)  which  project  onto  both  the  rotational  and  gravatational  modes  in  such  a  way  that 
transient  gravity  waves  are  minimized. 

In  Chapter  3  we  considered  the  formation  of  inner  bands,  which  we  defined  as  bands 
that  typically  form  less  than  500  km  from  the  center  of  the  vortex  and  owe  their  existence 
to  asymmetries  in  the  PV  field.  We  demonstrated  that  the  concept  of  PV  wave  break¬ 
ing,  which  has  been  found  to  be  useful  when  examining  stratospheric  data  (McIntyre  and 
Palmer,  1984;  Juckes  and  McIntyre,  1987;  Polvani  and  Plumb,  1992),  also  appears  to  be 
useful  for  explaining  the  formation  of  inner  hurricane  bands.  Simulations  of  this  process 
were  performed  by  assuming  convective  heating  had  generated  an  asymmetric  PV  field 
and  letting  the  model  evolve  unforced.  From  these  experiments,  we  found  that  the  PV 
contours  underwent  a  rapid  and  irreversible  deformation  which  acted  to  eject  filaments 
of  high  PV  air  from  the  center  of  the  vortex  to  large  radii,  thus  forming  bands.  Addi¬ 
tionally,  the  PV  field  also  underwent  a  continual  axisymmetrization  process  which  acted 
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to  restore  the  PV  field  to  a  more  circular  shape.  These  results  corresponded  well  to  the 
nondivergent  barotropic  results  found  by  Melander  et  a 1.  (1987b).  The  effect  of  the  wave 
breaking  and  axisymmetrization  process  was  to  cause  both  an  outward  shift  in  the  radius 
of  maximum  winds  and  an  increase  in  the  maximum  wind  speed.  Outside  the  radius  of 
maximum  winds,  but  inside  the  spiral  bands,  the  effect  of  symmetrization  was  to  decrease 
the  tangentially  averaged  wind  speed.  Near  the  banded  regions,  a  slight  increase  in  the 
tangentially  averaged  wind  speed  was  noticed. 

To  better  understand  the  wave  breaking  process,  passive  tracer  integrations  and  wave 
activity/wave  activity  flux  diagnostics  were  performed.  The  passive  tracer  experiments 
showed  that  the  PV  field  behaved  in  a  nearly  identical  fashion  to  a  passive  tracer,  even 
though  it  is  the  inversion  of  PV  which  determines  the  wind  and  mass  fields  for  balanced 
systems.  We  also  used  the  passive  tracer  technique  to  perform  Lagrangian  studies  of  the 
wave  breaking  process,  by  combining  two  passive  tracers,  one  which  initially  varied  only 
in  y  and  one  which  initially  varied  only  in  x.  Since  the  intersection  of  any  two  contours 
uniquely  defined  a  particle,  we  could  trace  the  trajectory  of  a  particle  simply  by  following 
the  path  of  the  proper  contour  intersection.  This  study  allowed  us  to  verify  that  the 
particles  which  appeared  in  the  banded  regions,  originated  near  the  lobes  of  the  elliptical 
shaped  PV  patterns. 

The  second  tool  used  to  study  the  wave  breaking  process  was  wave  activity  and  wave 
activity  flux  diagnostics.  The  exact  form  of  wave  activity,  as  derived  by  Haynes  (1988), 
was  used.  Wave  activity  provides  a  conservative  measure  of  wave  amplitude,  and  confirms 
that  banded  regions  are  also  regions  of  large  wave  amplitude.  The  implications  of  this, 
however,  are  unclear  since  the  form  of  wave  activity  used  is  incomplete  in  the  sense  that 
its  flux  cannot  be  directly  related  to  changes  of  the  mean  wind  (as  with  quasi-geostrophic 
theory).  We  therefore  cannot  accurately  determine  how  the  waves  and  mean  flows  will 
interact.  The  wave  activity  experiments  done  in  this  study  were  of  the  simplest  form; 
that  is,  there  was  no  external  forcing  involved  and  the  basic  state  was  both  axisymmetric 
and  constant  in  time.  It  is  expected  that  the  use  of  wave  activity  with  a  more  physically 
complete  model  would  prove  to  be  more  useful.  In  the  case  of  external  forcings,  the  wave 
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activity  would  not  be  conserved  in  diabatically  or  frictionally  forced  regions.  This  would 
allow  one  to  view  the  effects  of  external  forcings  on  wave  amplitude.  The  wave  activity 
studies  considered  here  were  merely  a  step  towards  this  goal. 

The  formation  of  bands  by  vortex  merger  was  also  studied.  Five  model  simulations 
were  presented  to  examine  the  banding  process.  The  first  four  simulations  considered  the 
merger  of  a  symmetric  vortex  with  pre-existing  PV  regions  of  different  sizes  and  ii_  jnsities. 
In  all  cases  spiral  bands  were  readily  formed.  In  the  fifth  simulation  we  considered  a 
symmetric  vortex  which  was  located  near  a  mass  sink  (PV  source).  The  mass  sink  was 
intended  to  represent  a  convective  heat  source  in  nature.  The  mass  sink  (heat  source)  was 
turned  on  in  a  gradual  manner  to  avoid  the  generation  of  transient  gravity  waves,  thus 
ensuring  the  model  evolved  along  the  diabatic  slow  manifold.  As  the  heating  intensified, 
PV  anomalies  were  induced  and  the  initial  vortex  propagated  towards  the  heat  source. 
The  effect  of  the  heat  source  was  to  generate  a  region  of  high  PV  air  which  was  advected 
around  and  into  the  vortex.  As  the  vortex  moved  closer  to  the  heat  source,  vortex  merger 
occurred  and  a  spiral  band  formed.  Unlike  the  unforced  experiments,  the  bands  associated 
with  the  forced  experiments  appeared  to  exhibit  deviations  from  axisymmetry  in  the  height 
fields.  These  results  relate  more  closely  to  observed  studies  (e.g.,  Ryan  et  a 1.,  1992;  Powell, 
1990)  and  suggest  that  hurricane  bands  evolve  along  the  diabatic  slow  manifold  (or  more 
precisely,  the  diabatic  frictional  slow  manifold)  in  nature. 

Results  from  Chapter  3  suggest  that  more  work  is  needed  to  fully  understand  hurri¬ 
cane  band  formation  by  wave  breaking  and  vortex  merger  processes.  At  the  observational 
level,  PV  maps  of  tropical  cyclones  are  desperately  needed.  This  is  not  an  easy  task,  how¬ 
ever,  since  the  calculation  of  Rossby  Ertel  PV  (which  can  be  approximated  by  —  g( (86/ dp)) 
involves  the  calculation  of  both  horizontal  and  vertical  derivatives.  Any  errors  in  the  wind, 
are  multiplied  when  their  derivatives  are  calculated.  In  addition,  obtaining  a  sufficiently 
dense  data  set  to  resolve  spiral  bands  (~30  km  wide)  is  difficult  for  obvious  reasons.  The 
use  of  aircraft  Doppler  radar  in  recent  years,  however,  has  shown  promise  by  provid¬ 
ing  a  relatively  accurate  and  dense  source  of  wind  data.  Work  to  compute  PV  fields  in 
hurricanes  is  currently  being  attempted  at  the  Hurricane  Research  Division  of  NOAA. 
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Aside  from  observational  studies,  more  work  in  the  area  of  modeling  is  also  needed. 
The  results  presented  in  Chapter  5  were  for  the  simplest  case  possible.  Although  this 
allowed  us  to  isolate  the  fundamental  dynamics,  other  factors  related  to  tropical  cyclones 
may  play  a  key  role.  The  foremost  factor  is  friction.  Since  the  model  used  in  this  paper 
was  a  shallow  water  model,  frictional  forces  could  not  be  represented  in  any  true  physical 
sense.  That  is,  we  could  add  a  surface  friction  term  in  the  form  of  a  bulk  aerodynamic 
formula,  but  its  effects  would  be  to  simply  damp  the  winds  with  time  -  no  feedback  into  a 
cumulus  parameterization  scheme  is  possible.  This  suggests  that  a  layered  normal  mode 
model  would  also  be  useful  since  frictional  effects  could  then  be  fed.  back  into  the  heating 
fields.  An  incompressible,  three  layer  (fixed  boundary  layer,  inflow  region,  and  outflow 
region)  normal  mode  model  of  this  type  was  accomplished  by  DeMaria  and  Schubert 
(1984)  and  was  used  to  study  tropical  cyclone  motion.  A  similar,  but  more  useful  normal 
mode  model,  which  uses  three  isentropic  layers  was  developed  by  DeMaria  and  Pickle 
(1988)  and  was  used  by  Shapiro  (1992)  to  investigate  the  effect  of  asymmetries  and  the 
/3-effect  on  hurricane  motion.  Preliminary  results  towards  this  goal  have  been  attempted 
in  this  study  through  use  of  a  three  layer  extension  of  the  present  model.  The  three  layers 
represented  the  boundary  layer,  inflow  (cyclonic)  layer,  and  outflow  (anticyclonic)  layer 
of  a  tropical  cyclone.  In  contrast  to  DeMaria  and  Schubert  (1984),  the  boundary  layer 
was  of  variable  depth.  When  the  model  was  run  unforced  with  an  initial  elliptical  shaped 
PV  pattern,  results  were  similar  to  those  obtained  for  the  shallow  water  case.  Additional 
work  in  this  area  is  still  needed  to  add  a  workable  cumulus  parameterization  scheme. 

In  Chapter  4,  we  considered  the  formation  of  outer  spiral  bands.  These  bands  were 
defined  as  those  which  typically  form  500  km  from  the  vortex  center,  and  owe  their  exis¬ 
tence  to  the  breakdown  of  the  ITCZ.  Several  general  stability  theorems,  both  linear  and 
nonlinear,  were  reviewed.  These  theorems  all  indicated  that  a  reversal  in  the  PV  gradient 
is  a  necessary  condition  for  barotropic  instability.  Although  the  general  stability  theo¬ 
rems  provided  conditions  for  instability,  they  offered  no  information  concerning  growth 
rates  or  unstable  patterns.  In  order  to  predict  such  features,  we  utilized  the  less  general 
technique  of  linear  normal  mode  stability  analysis.  This  analysis  technique  is  less  general 
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because  it  assumes  the  instabilities  grow  with  spatial  structure  e,kx.  Results  from  the 
normal  mode  stability  analysis  showed  that  as  the  width  of  a  PV  anomaly  decreased,  the 
growth  rates  of  the  instabilities  increased  and  the  wavenumber  of  the  most  unstable  mode 
also  increased.  To  test  the  results  of  the  normal  mode  stability  analysis,  numerous  nu¬ 
merical  model  integrations  were  performed  by  initializing  the  model  with  various  shapes 
and  sizes  of  PV  strips.  These  strips  were  continuous  in  x  but  confined  in  y  except  for 
the  final  experiment.  In  the  final  experiment,  the  initial  PV  field  was  confined  in  both  x 
and  y.  Results  from  the  numerical  simulations  displayed  similar  trends  to  those  predicted 
by  linear  theory.  That  is,  instabilities  grew  faster  for  the  thinner  strips  and  showed  a 
tendency  for  higher  wavenumbers.  In  addition,  these  results  bore  remarkable  resemblance 
to  the  natural  breakdown  of  the  IT^Z  over  the  eastern  Pacific  Ocean.  One  implication  of 
these  results  is  that  the  analysis  of  easterly  wave  formation  in  the  North  African  region 
performed  by  Burpee  (1972)  is  not  unique  to  that  region.  As  Schubert  et  al.  (1991) 
demonstrated,  ITCZ  convection  alone  can  lead  to  PV  anomalies.  These  PV  anomalies, 
in  turn,  satisfy  the  necessary  condition  for  unstable  flows,  resulting  in  breakdowns  in  the 
ITCZ  similar  to  those  presented  in  Chapter  4. 

The  resulting  patterns  from  initially  similar  PV  strips  appears  to  be  highly  sensitive 
to  the  initial  condition.  This  sensitivity  was  tested  by  making  three  separate  model  inte¬ 
grations  starting  from  identically  specified  wind  and  mass  fields.  In  each  case,  the  fields 
were  perturbed  using  different  sequences  of  random  numbers  (but  of  equal  magnitude). 
In  all  three  cases  the  ITCZ  broke  down  and  pools  of  high  PV  air  were  formed.  The  dif¬ 
ference,  however,  was  in  the  number  of  pools  which  formed.  The  three  runs  resulted  in 
one  wavenumber  three  pattern  and  two  wavenumber  two  patterns.  This  suggests  that 
predicting  the  locations  of  incipient  hurricanes  based  on  ITCZ  patterns  would  be  difficult 
because  of  the  inherent  nonlinearity  of  the  problem. 

One  possible  extension  of  the  previous  study  is  to  include  spherical  effects.  Some 
work  in  this  area  has  been  performed  by  Dritschel  (1992)  using  a  nondivergent  barotropic 
model  on  a  sphere.  He  concluded  that  barotropic  flows  on  a  sphere  have  an  even  more 
pronounced  tendency  to  produce  small  short  lived  vortices,  especially  in  equatorial  and 
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mid-latitude  regions.  Work  is  currently  being  done  at  C.S.U.  to  extend  Dritschel’s  work 
to  the  shallow  water  system  using  the  recently  developed  NCAR  global  shallow  water 
spectral  model  (Hack  and  Jakob,  1992). 

In  Chapter  5,  we  investigated  the  stability  of  spiral  bands.  This  investigation  followed 
logically  from  Chapters  3  and  4,  the  reason  being  that  in  Chapter  3  we  saw  how  the  banding 
process  created  strips  of  relatively  high  PV  (and  therefore  reversals  in  the  PV  gradient), 
and  in  Chapter  4  we  saw  that  regions  of  PV  gradient  reversals  satisfied  the  necessary 
condition  for  barotropic  instability. 

The  stability  analysis  was  done  both  analytically  and  numerically.  The  analytic 
study  again  made  use  of  nonlinear  normal  mode  analysis.  An  eigenvalue/eigenvector 
problem  was  formed  by  considering  the  boundary  displacement  of  multiple  concentric 
annular  regions  of  uniform  vorticity.  To  simulate  an  isolated  spiral  band,  we  considered 
the  case  of  a  single  annular  region  of  PV  located  at  an  arbitrary  distance  from  the  center 
of  the  domain.  The  stability  analysis  produced  similar  results  to  that  of  Chapter  4.  As  the 
band  was  made  narrower,  both  the  growth  rates  and  wavenumbers  of  the  most  unstable 
mode  increased.  These  results  were  also  verified  through  numerical  integrations. 

The  above  mentioned  numerical  simulations  lead  to  a  somewhat  unexpected  finding. 
As  the  annular  regions  of  PV  broke  down  in  the  above  experiments,  they  began  to  form 
polygonal  shapes.  Similar  polygonal  shapes  have  also  been  found  in  the  eyewalls  of  intense 
hurricanes  (Lewis  and  Hawkins,  1982;  Muramatsu,  1986).  If  we  approximate  hurricane 
eyewalls  by  annular  regions  of  high  PV  air,  then  it  seems  reasonable  to  assume  that 
barotropic  instabilities  should  also  be  observed,  especially  in  more  intense  hurricanes. 
As  a  result,  these  barotropic  instabilities  could  act  to  form  polygonal  shapes  in  hurricane 
eyewalls.  As  further  observational  evidence  of  eyewall  instabilities,  Black  and  Marks  (1991) 
noted  the  presence  of  a  meso- vortex  in  Hurricane  Hugo  (1989).  This  may  be  due  to  the  PV 
pooling  process  which  occurs  as  instabilities  grow  Numerical  simulations  demonstrated 
well  the  breakdown  of  an  annular  PV  region  and  subsequent  formation  of  smaller  scale 
circulations.  On  a  much  finer  scale,  the  instabilities  of  annular  high  PV  regions  may 
also  have  some  relevance  to  the  observed  phenomenon  of  tornado  suction  spots.  These 
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suction  spots  are  secondary  vortices  which  rotate  around  the  parent  vortex  (Fujita  et 
a i.,  1972).  Linear  normal  mode  analysis  of  this  phenomena  has  been  performed  by  Snow 
(1978),  and  bear  some  resemblence  to  our  present  results.  Although  these  two  phenomena 
may  be  closely  related,  tornado  suction  spots  are  obviously  a  nonhydrostatic  process  and 
therefore  not  representable  by  a  shallow  water  model  such  as  the  one  used  in  the  present 
study. 

In  the  final  section  of  Chapter  5,  we  offered  an  explanation  as  to  why  spiral  bands  are 
not  observed  to  undergo  unstable  growth  even  though  they  satisfy  the  necessary  condition 
for  instability.  The  explanation  appears  to  be  related  to  the  adverse  shear  generated  by 
the  main  vortex.  Physically,  the  adverse  shear  of  the  main  vortex  opposes  the  circulation 
which  the  PV  anomaly  tries  to  generate.  Results  which  demonstrate  this  effect  in  a 
two-dimensional  fluid  have  been  presented  both  numerically  and  analytically  by  Dritschel 
(1989).  In  the  present  study,  we  also  used  both  linear  normal  mode  analysis  and  numerical 
simulations  to  study  the  problem.  The  linear  analysis  examined  the  stability  of  a  50  km 
wide  annular  region  of  PV  centered  at  225  km  from  the  cylindrical  coordinate  center. 
At  the  coordinate  center  was  a  uniform  vorticity  region  of  50  km  radius.  The  vorticity 
of  the  center  region  was  allowed  to  vary,  thus  providing  different  intensities  of  adverse 
shear.  Results  showed  that  as  the  adverse  shear  increased,  the  growth  rate  of  the  most 
unstable  mode  decreased,  while  the  wavenumber  increased.  A  similar  case  to  that  used 
for  the  linear  analysis  was  simulated  using  the  model.  In  the  case  of  the  model  simulation, 
the  vorticity  in  the  main  vortex  was  twice  that  of  the  band.  Results  from  the  simulation 
confirmed  those  suggested  by  the  linear  analysis.  The  complete  breakdown  of  the  annular 
region  took  approximately  three  times  longer  and  the  wavenumber  of  the  most  unstable 
mode  doubled. 

A  problem  that  persisted  throughout  the  Chapter  5  experiments  was  the  tendency 
for  the  annular  PV  regions  to  be  perturbed  in  a  wavenumber  four  pattern.  This  tendency 
occurred  because  of  doubly  periodic  boundary  conditions.  To  overcome  this  problem,  it 
may  be  useful  to  integrate  the  cylindrical  coordinate  model  presented  in  Appendix  A. 
This  woula  allow  random  perturbations  in  the  initial  fields  to  be  used  for  determining  the 
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true  most  unstable  mode.  Two  drawbacks,  however,  are  that  the  basis  functions  of  the 
cylindrical  model  have  poor  convergence  properties,  and  no  fast  transform  exists  in  the 
radial  direction.  For  simple  studies,  such  as  those  suggested  above,  these  problems  may 
be  surmountable. 
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Appendix  A 


THE  NORMAL  MODE  MODEL  IN  CYLINDRICAL  COORDINATES 


In  this  chapter  we  present  the  normal  mode  transform  in  cylindrical  coordinate.  Al¬ 
though  this  version  of  the  model  was  not  used  in  the  numerical  integrations  it  is  none  the 
less  useful  to  see  the  similarities  to  the  Cartesian  coordinate  model.  It  is  of  interest  to 
note  that  all  normal  mode  models  reduce  to  the  oscillation  equation  (2.18)  regardless  of 
geometry  or  vertical  structure. 

A.l  The  Governing  Equations  in  Cylindrical  Coordinates 


The  shallow  water  equations  in  cylindrical  coordinates  can  be  written  in  rotational 
form  as 


~-f(v  +  A)  +  ^  +  K)  =  F 

(A.l) 

f +'<t,+B>  +  ;J^+*>  =  G 

(A.2) 

d<f>  ,  2  fdr(u  +  C)  ,  d(v  +  D)\  „ 

dt  +  V  rdr  rd<p  )  ’ 

(A.3) 

where  u  is  now  the  radial  component  of  velocity,  v  is  the  tangential  component,  and 


n  _  ^  d(rv)  du 

/’  ^  rdr  rdtp 


K  =  ±(u2  +  v2) 


C 


n  _ 

H'  H' 


and  F,  G,  and  5  are  the  frictional  and  diabatic  effects.  Using  these  equations  we  can  now 
complete  the  tangential  transform 
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A.2  The  Tangential  Transform  and  Lateral  Boundary  Condition 

If  we  define  the  tangential  finite  Fourier  transform  of  u(r,  <p,  t)  as  un(r,  t),  the  tangen- 


tial  transform  pair  can  be  written  as 

un(r,  t)  =  u(r,  ip,  t)e~>nifidip 

(A. 4) 

u(r,<p,t)=  52  «n (r,*)e,nv 

(A.5) 

n=  —  oo 


where  n  represents  the  wavenumber  in  the  tangential  direction.  By  using  similar  definitions 
for  the  tangential  win  1  component  and  geopotentia) ,  VA.1)-(A.3)  can  be  transformed  to 
obtain 


^  -  f(vn  +  An)  +  !-(<*>„  +  Kn)  =  F„ 

(A. 6) 

%  +  /K  +  Bn)  +  -(*„  +  Kn)  =  Gn 
at  r 

(A. 7) 

d<t>n 

dt 

+  ~(vn  +  A»)j  =  Sn , 

(A.8) 

Having  completed  the  tangential  transform,  we  are  now  far  enough  into  spectral  space 
co  define  our  lateral  boundary  condition.  Since  their  is  no  natural  physical  boundary 
condition,  we  will  choose  the  simplest  computational  boundary  condition.  We  consider 
the  perturbation  geopotential,  <f>,  to  vanish  at  r  =  a  for  all  time  (where  a  is  the  radial 
boundary).  If  we  also  assume  that  the  nonlinear  and  forcing  terms  vanish  for  large  r, 
then  from  (A. 3)  the  divergence  must  also  vanish  at  the  boundary.  Thus,  our  boundary 
condition  is  nearly  equivalent  to  a  zero  di  *rgence  boundary  condition. 

Before  considering  the  problem  of  radially  transforming  our  governing  equations,  it 
will  again  prove  convenient  to  express  our  tangentially  transformed  equations  in  vector 
form.  We  can  accomplish  this  as  we  did  in  chapter  three  by  defining  W„  —  [u„,  vn,  <pn/c)T . 
Again  these  vectors  have  been  chosen  to  ensure  consistency  in  units.  Using  these  definitions 
we  can  write 

dWn 


dt 


+  £nW„  =  Fn  +  N„, 


(A. 9) 
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where  the  nonlinear  terms  Nn  are  given  as 

N„  =  [m.  -  £k„,  -fBn  -  ±K.,  -c  +  ^a,)]T  (A.10) 

and  the  Forcing  terms,  Fn,  are  given  as 

F „  =  [F„,  Gn,  Sn/c]T.  (A. 11) 

The  matrix  operator  is  defined  as 

1  0  “/  ch  N 

Cn  =  /  0  cf  .  (A. 12) 

c?  0 

Using  these  definitions,  we  are  now  prepared  to  complete  the  radial  transform  of  (A. 9). 

A.3  The  Radial  Transform 

We  begin  the  radial  transform  by  defining  the  complex  inner  product  to  be 

2  [a 

(X,  y )  =  ^2  JQ  (*0«o  +  *iy[  +  x2V2)  rrfr-  (A. 13) 

where  x  and  y  are  three  component  complex  vectors  which  satisfy  the  boundary  condition 
and  the  asterisk  denotes  the  complex  conjugate.  As  discussed  in  Appendix  B,  £„  is  skew- 
Hermitian  with  respect  to  the  inner  product  (A. 13),  meaning  £j,  =  — £„.  Since  it  is  well 
known  that  the  eigenvalues  of  skew-Hermitian  operators  are  pure  imaginary,  we  can  write 

£? jK„pg  =  ilSnpqlCnpq,  (A. 14) 

where  i/npq  are  the  eigenvalues  and  Knp9  are  the  corresponding  vector  eigenfunctions  (or 
transformation  kernels).  These  subscripts  will  become  clearly  defined  shortly.  One  cam 
easily  solve  for  these  eigenvalues  and  eigenfunctions  by  substituting  our  definition  of  £n 
into  (A.  14)  and  reducing  the  system  of  three  equations  down  to  one  equation  in  one 
unknown.  The  equation  obtained  by  solving  for  the  third  component  of  Knp9  is  the 
order  n  Bessel  equation  and  the  corresponding  solution  is  chosen  to  be  the  order  n  Bessel 
function  of  the  first  kind,  Jn(k„pr ),  where  k%p  =  (v £  —  f2)/c It  should  be  noted  that 
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we  chose  Bessel  functions  of  the  first  kind  to  ensure  our  solutions  remain  bounded  at  the 


origin.  The  other  two  components  are  then  obtained  through  direct  substitution. 

One  other  possible  solution  for  K npq  can  also  be  obtained.  This  corresponds  the  the 
case  when  vnpq  =  0.  We  will  again  refer  to  this  as  the  geostrophic  case  since  the  equations 
which  are  obtained  by  substituting  vnpq  =  0  into  our  definition  of  Cn  are  simply  the 
geostrophic  balance  equations.  It  is  important  to  note  that  any  function  can  be  chosen 
for  the  third  component  of  the  kernel  and  it  will  satisfy  the  geostrophic  relation.  We 
again  choose  the  Bessel  function  of  the  first  kind  to  remain  consistent  with  our  previous 
solutions.  This  case  will  be  denoted  with  the  subscript  q  =  0. 

Finally,  we  can  write  the  transformation  kernel  as 


-'j?Jr>(knpr) 

cdJn(k.rr) 


g  =  0 


\  Jn(knpr)  ) 

Knpq(r)  =  Anpq  « 

(  \ 

?fiZMk.pr) ,  =  1,2 

.  V  ^fT-MknpT)  J 

where  Anpq  is  the  normalization  factor  given  by 

Anpq  =  /  {(^M  +  /2  +  c2fc2p)J2+1(fcnPa)}  2  . 


(A.15) 


(A. 16) 


Jn  is  the  order  n  Bessel  function  of  the  first  kind  and  knp(p  —  1,2,3,...)  are  the  discrete 
wave  numbers  satisfying  the  boundary  condition 


Jn(,knpa)  —  0, 


(A. 17) 


0  9  =  0  (geostrophic  modes) 

i/nM  =  '  ~(/2  +  c2^)!  9  =  1  (gravity  -  inertia  mode)  >.  (A. 18) 

.  (/2  +  ^np)*  9  =  2  (gravity  -  inertia  mode)  , 

Equation  (A.  17)  is  derived  by  setting  the  third  component  of  Knpq  equal  to  zero  at  the 
boundary  in  order  to  satisfy  the  boundary  condition  that  <f>  must  vanish  at  r  =  a.  Since 
/,  c,  and  k  are  all  non-zero,  the  only  possible  solution  is  for  Jn(knpa)  to  be  zero  for  both 
the  geostrophic  and  gravity-inertia  modes. 
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It  should  now  be  clear  just  exactly  what  the  subscripts  n,  p,  and  q  represent.  The  n 
obviously  refers  to  the  tangential  mode,  the  p  indexes  the  radial  modes,  and  the  q  gives 
the  mode  type  corresponding  to  the  three  possible  solutions  for  the  frequency. 

Before  we  discuss  the  orthonormality  of  our  basis  functions,  we  must  first  consider 
how  they  behave  at  the  origin.  We  know  that  our  kernel  must  be  bounded  at  the  origin 
(because  we  chose  Bessel  functions  of  the  first  kind),  but  we  do  not  know  what  values  the 
kernels  take  at  the  origin.  In  order  to  examine  this,  we  use  the  limiting  form  of  Bessel 
functions  of  the  first  kind  for  small  arguments  (Abramowitz  and  Stegun,  eqn.  9.1.7),  that 


where  Wnpq(t)  is  a  complex,  scalar  function  of  time  having  units  of  ms-1.  Equation 
(A.23)  is  the  transform  from  radial  physical  space  to  spectral  space  and  (A. 22)  is  the 
inverse  transform. 
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To  derive  the  spectral  equations  (i.e.  equations  involving  Wnpq),  we  take  the  inner 
product  of  (A.9)  with  Knpq(r)  to  obtain 

^(Wn,Knw)  +  (£„Wn,Knw)  =  (Fn,K„M)  +  (Nn.Knp,).  (A. 24) 

Since  £),  =  -£„  with  respect  to  the  inner  product,  (A. 13)  we  can  use  (A.12)  to 
rewrite  the  second  term  in  (A.24)  as 

(Cn  w„,  Knpq)  =  -(W„,  CnKnpq)  =  -ti/„„(Wn,  Knpq).  (A.25) 

Using  (A.24)  we  see  that  (A.25)  can  be  written  as 

dW 

21  +  i»npqWnpq  =  Fnpq  +  Nnpq.  (A. 26) 

Equation  (A. 26)  completes  our  transformation  from  our  original  governing  equations  in 
physical  space  to  normal  mode  spectral  space.  It  is  of  interest  to  note  that  any  normal 
model  spectral  model  (regardless  of  the  equations,  geometry,  etc.)  leads  to  an  equation  of 
this  form. 

A.4  Computational  Procedures 

This  section  discusses  some  the  complications  which  arise  when  integrating  the  cylin¬ 
drical  coordinate  model.  Included  are  discussions  of  the  horizontal  transforms  and  the 
calculation  of  derivatives  and  nonlinear  terms. 

A.4.1  The  Horizontal  Transforms 

The  tangential  transform  can  be  easily  accomplished  with  any  real  FFT  routine  in 
a  similar  fashion  to  the  transforms  discussed  in  chapter  three.  Since  FFTs  are  well  doc¬ 
umented  in  academic  texts  (e.g.  Brigham,  1988),  the  tangential  transform  will  not  be 
discussed  further  here.  The  radial  transform,  however,  does  require  some  discussion. 

Unlike  the  Cartesian  coordinate  model,  our  transformation  kernel  (A. 15)  does  not 
consist  of  periodic  functions.  Because  of  this,  there  is  no  hope  of  computing  the  radial 
transform  (A.  11)  exactly  with  a  finite  number  of  points.  Aliasing  error  is  therefore  un¬ 
avoidable. 
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Another  difficulty  with  the  radial  transform  is  that  there  exists  no  fast  transform  for 
Bessel  functions.  We  are  therefore  forced  to  evaluate  the  integral  in  (A.  13)  using  the  full 
number  of  operations  required  by  the  chosen  quadrature  method. 

In  the  present  study,  we  computed  the  integral  in  (A. 13)  using  3P+1  points  using  both 
trapezoidal  and  Simpson’s  quadrature,  where  P  represents  the  subscript  of  the  highest 
radial  wave  number.  We  chose  3 P  + 1  points  in  order  to  be  consistent  with  the  Cartesian 
coordinate  model.  Results  from  linear  model  integrations  indicated  that  with  this  number 
of  points,  Simpson’s  quadrature  produced  the  best  results. 

Finally,  it  should  be  noted  that  Bessel  series  have  a  relatively  slow  rate  of  convergence 
which  is  on  the  order  of  1/N  where  N  is  the  number  of  terms  in  the  series  (Gottlieb  and 
Orzag,  1977).  Fourier  series  on  the  other  hand  converge  on  the  order  of  1/N2  (Gottlieb 
and  Orzag,  1977).  This  implies  that  more  waves  would  be  required  for  the  Bessel  series 
to  have  the  same  accuracy  as  the  Fourier  series  for  the  same  function. 

A.4.2  Calculation  of  Derivatives  and  Nonlinear  Terms 

It  can  easily  be  shown  from  (A.22)  that  r-derivatives  can  be  calculated  by  summing 
the  spectral  coefficients  Wnpq  with  the  r-derivatives  of  the  transformation  kernel.  Since 
derivatives  of  Bessel  functions  are  known  exactly  (see  B.9),  derivatives  of  the  transforma¬ 
tion  kernel  axe  known  exactly.  This  indicates  that  derivatives  of  the  dependent  variables 
can  be  computed  exactly  to  within  the  truncation  limits  of  the  model.  One  drawback 
however  is  that  the  transformation  kernels  are  not  eigenfunctions  of  the  d/dr  operator. 
Thus,  derivatives  of  the  transformation  kernels  must  be  calculated  at  each  time  step  or  the 
derivative  of  the  transformation  kernel  must  be  stored  in  an  array.  The  choice  of  course 
depends  on  the  amount  of  available  computer  storage.  To  demonstrate  the  amount  of 
storage  required  for  transformation  kernels,  consider  the  case  where  we  use  128  points  in 
the  r  direction  to  evaluate  the  transformation  kernel,  if  we  keep  42  waves  in  both  the  radial 
and  tangential  directions  then  roughly  5.4  megabytes  of  storage  are  required  (assuming  a 
complex  variable  is  stored  as  two  four  byte  words).  If  we  choose  to  double  this  resolution, 
the  amount  of  required  storage  increases  by  a  factor  of  eight  to  43.3  megabytes.  One  can 
easily  see  that  as  the  resolution  improves,  the  model  quickly  becomes  cost  prohibitive. 
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Nonlinear  terms  would  require  enormous  computational  expense  in  cylindrical  coor¬ 
dinates.  The  reason  being  that  not  only  would  a  transform  to  and  from  normal  mode 
spectral  space  be  required  at  each  time  step,  but  derivatives  of  the  transformation  kernel 
would  need  to  be  stored  or  calculated. 


Appendix  B 


THE  SKEW-HERMITIAN  PROPERTIES  OF  THE  LINEAR  OPERATORS 


In  this  appendix  we  will  prove  the  skew-Hermitian  property  of  the  matrix  operators 
(2.5)  and  (A.12)  with  respect  to  the  inner  products  (2.8)  and  (A. 13)  respectively.  Our 
first  step  is  to  define  what  is  meant  by  skew-Hermitian.  A  linear  operator  is  said  to  be 
skew-Hermitian  with  respect  to  an  inner  product  if  its  adjoint  is  the  same  operator  with 
the  opposite  sign.  We  define  the  adjoint  operator,  C\  to  be  the  operator  which  satisfies 

(£p,q)  =  (p,£tq),  (B.l) 

where  p  and  q  in  (B.l)  represent  complex  vectors  which  satisfy  the  conditions  of  the  inner 
product.  Therefore,  if  an  operator  is  skew-Hermitian  with  respect  to  an  inner  product, 
then 

(£p,q)  =  (p,-£q).  (B.2) 

B.l  The  Cartesian  Coordinate  System 


Consider  p  and  q  to  be  the  complex,  three  component  vectors  which  are  periodic  in 
both  x  and  y  with  period  Lx  and  Lv  respectively.  Using  (2.5)  and  (2.8)  we  can  write 


(c 


,  ,  „<9po  ,  dpi' 

+  l  a7  +  Cdy 


(B.3) 


q^dxdy 

If  (B.3)  is  integrated  by  parts  and  terms  are  rearranged,  it  can  be  rewritten  as 

<£Pi  q) = ihr,  r  t  (/,r  “ '  & ) + 1 n  (-/,J  -  CW + 

n  ("e^ ' c^) +  c  (s<P2,J)  +  + slpo,;) + |;(pi,5)) } dxdy  <B-4) 
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Because  both  p  and  q  periodic  in  x  and  y,  the  last  term  in  (B.4)  must  be  identically  zero. 
By  factoring  the  complex  conjugate  operator  from  terms  in  parentheses,  we  see  that  (B.4) 
reduces  to 

(£p,q)  =  (p,-£q).  (B.5) 

Since  =  -£,  we  conclude  that  £  is  skew-Hermitian  with  respect  to  (2.8). 

B.2  The  Cylindrical  Coordinate  System 

Now  consider  p  and  q  to  be  the  complex,  three  component  vectors  whose  third  com¬ 
ponents  vanish  at  r  =  a.  Using  (A.12)  and  (A. 13)  we  can  write 

(£„p,  q)  =  ^  jT  |  (~/pi  +  )  qo  +  (jpo  +  c^p2)  q*i 

+  (" 't£t  +  cTPl) q 5 } rdr ■  (B6) 

If  we  integrate  by  parts  the  terms  involving  derivatives  with  respect  to  r  and  rearrange 
terms  we  get 

O 

(£„p,q)  =  ^  yjpo  ( fq\  ~  c^)  +Pl  (-/« 0  + 

0 

+P2  }  rdr  +  [flP072  +  aP2?olr=a  ■  (B-7) 

We  know  from  our  boundary  condition,  however,  that  p2  and  q%  are  zero  at  r  =  a  therefore 
the  boundary  terms  in  (B.7)  must  vanish.  In  addition,  if  we  factor  out  the  complex 
conjugate  operator  from  the  terms  in  parentheses,  (B.7)  can  be  written  as 

(£„p,  q)  =  (p,  -£nq).  (B.8) 

Since  C\  =  —  £„,  we  again  conclude  that  £„  is  skew-Hermitian  with  respect  to  the  inner 
product  (A. 13). 


Appendix  C 


THE  ORTHONORMALITY  OF  THE  BASIS  FUNCTIONS 

In  this  appendix  we  prove  the  statement  of  orthonormality  of  our  basis  functions 
(2.13)  and  (A.15)  with  respect  to  the  inner  products  (2.8)  and  A.13)  respectively.  The 
proof  will  be  done  first  in  Cartesian  coordinates  and  then  in  cylindrical  coordinates. 

C.l  Orthonormality  of  Cartesian  Coordinate  Basis  Functions 

The  proof  of  orthonormality  in  Cartesian  coordinates  is  straightforward.  Since  the 
x  and  y  dependency  of  the  basis  functions  is  periodic,  it  is  easily  verified  that  if  k  ^  k' 
and/or  l  ^  l1,  the  integral  in  (2.8)  will  be  identically  zero.  If  the  k  =  k'  and  l  ~  I'  it  is  easily 
verified  that  the  x  and  y  dependency  will  integrate  to  exactly  one,  leaving  only  a  vector 
multiplication  left  to  evaluate.  Because  of  this,  if  we  wish  to  show  that  the  eigenfunctions 
are  orthonormal,  we  need  only  to  show  that  the  x  and  y  independent  parts  of  (2.13)  are 
orthonormal.  One  can  easily  prove  this  to  be  true  by  computing  the  vector  products  of 
the  as  defined  in  (2.8)  for  the  x  and  y  independent  parts  of  (2.13)  for  each  mode  type,  q. 

C.2  Orthonormality  of  Cylindrical  Coordinate  Basis  Functions 

The  proof  of  orthonormality  in  cylindrical  coordinates  is  considerably  less  straight¬ 
forward.  Our  first  step  is  to  prove  the  orthogonality  of  Knpq  with  Knpy  where  p^p',q  = 
1  or  2,  and  <f  =  0,  1,  or  2.  We  do  this  by  adding  the  inner  product  of  CnKnpq  =  ivnj>qKnm 
with  Knj>Y  to  the  inner  product  of  Knpq  with  £nKnpy  =  wnp>q'Knpigi  to  obtain 

(^nKnpq,  Knp<q')  +  (Knpq,  £nKnp>qi)  — 


(iVnpqK-npqi  Knj/^)  +  (K  npq,  iVnpiq'  K  \p'q'  )  ■ 


(C.l) 
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Because  of  the  skew-Hermetian  property  of  £„  (as  discussed  in  appendix  A),  the  left  hand 
side  of  (C.l)  vanishes.  We  can  therefore  rewrite  (C.l)  as 


(^n pq  l'nptq' )  ,  K np'q' )  —  O' 


(C.2) 


Since  the  eigenvalues  for  the  gravity-inertia  modes  are  distinct,  the  corresponding  eigen¬ 
functions  must  be  orthogonal  to  satisfy  (C.2).  In  addition,  since  the  eigenvalues  for 
the  gravity-inertia  modes  are  non-zero,  the  inner  product  of  an  eigenfunction  from  the 
geostrophic  mode  set  must  be  orthogonal  to  an  eigenfunction  from  the  gravity-inertia 
mode  set  to  satisfy  (C.2). 

Now  we  wish  to  prove  the  orthogonality  of  the  degenerate  (geostrophic)  modes.  Con¬ 
sider  the  case  where  p  p'  and  q  —  q1  =  0.  The  inner  product  of  these  eigenfunctions 
is 

2  f 

(KnpOi  Knp<o)  =  |  ~f2^2  Jn(knp1')Jn(knj/r)~l~ 

j +  Jn(fcnpr)  Jn(fcn^r)  1  rtfr.  (C.3) 

If  we  integrate  the  second  term  on  the  right  hand  side  by  parts,  we  can  shift  the  r  derivative 
off  Jn{knpr)  and  rewrite  (C.3)  as 


(Knp0»  Knp»o)  —  AnpoAnp'O 


i  rl&Lj 

a2  JQ  |  /2r2 


(knpr)Jn(knpiT) 


~  pJn(knpr)rdr 


rdJn^T)  +  Jn{knpr)Jn(knprr)^  rdr+ 


T  ^ dJn(knf,r ) 

O  -p Jn  ( knpa) - ^ — 


(C.4) 


Equation  (C.4)  can  be  simplified  further  by  first  recognizing  that  the  boundary  term 
vanishes  because  of  (A. 17).  Secondly,  we  may  eliminate  the  derivative  in  the  second  term 
one  the  right  by  using  the  equation  that  Jn(knf/r)  satisfies  (i.e.  Bessel’s  equation)  divided 
by  r2.  This  equation  cam  be  written  as 

-.2 


^|rw:^r)]+(^_^)Jn(Vr)  =  0 


(C.5) 
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By  substituting  (C.5)  into  (C.4)  wo  obtain 


1  2  fa 

(KijpO?  Knp»o)  =  ■dnpO-dnp'Oyj  (^npO  4"  /  4"  cm^“np)  ^2  ^  Jn(knpr)Jn(knptr)rdr,  (C.6) 

where  i/npo  is  identically  zero  but  was  added  for  consistency  in  notation.  From  Abramowitz 
and  Stegun  (1964,  page  485,  eqn.  11.4.5),  the  integral  in  (C.6)  is  evaluated  as 


l  J.MW)*  - 1  [(! .  ^  jl{KtA) 


'  dJn(knrr)  \  21 


p=p' 


However,  this  can  be  simplified  to 

fa  (  0  p^p'  } 

/  Jn(fcnPr)7n(A:Bp,r)rdr={  .  (C.8) 

Jo  {  \J^{knpa)  p  =  p'  ) 

through  use  of  our  boundary  condition  (2.22)  and  the  derivative  relation  (Abramowitz 
and  Stegun,  1964,  page  361,  eqn.  9.1.30) 

-  kjn+1(kr).  (C-9) 

Finally,  by  substituting  (C.8)  into  (C.6)  while  using  the  (A. 16),  we  can  write 


f  0  p#p'  1 

(KnpOi  Kn^o)  =  {  f  • 

l  1  p=p'  J 


(C.10) 


Thus,  our  eigenfunctions  are  orthonormal  even  for  the  degenerate  (geostrophic)  case. 

Our  final  task  is  to  show  that  eigenfunctions  within  the  gravity-inertia  mode  sets  are 
orthonormal.  To  do  this,  consider  the  case  where  p  =  p'  and  q  =  q'  =  1  or  2.  Our  inner 


product  can  then  be  written  as 


v npq  f  dJn(knpr) ' 
f2k2  \  dr 


(Knp^Knp,)  -  AifM)a2  j  ^iT^n^V)  +  ^  ?  )  + 

+  jr (— ir^1)2 - +  ^.W)} -*■  (C.11) 
The  fifth  term  on  the  right  hand  side  of  (C.ll)  must  vanish  at  the  upper  boundary  because 
of  (A.  17)  and  also  at  the  lower  boundary  because  7„(0)  =  0  for  all  n  >  0.  We  can  also 
simplify  (C.ll)  as  we  did  for  the  geostrophic  modes  by  integrating  the  by  parts  the  terms 
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involving  derivatives  and  using  the  equation  which  Jn(knpr)  satisfies.  By  doing  this  and 
again  using  (A. 17),  we  can  rewrite  (C.ll)  as 

(Knpq,  Knpq)  ==  Anwy2  (^npq  /  d"  Cm^np)^2  *^n(^npr)r^r  (C.12) 

Using  (C.8)  and  (A.  16),  (C.12)  can  be  written  in  final  form  as 

(Knpq,  K„pq)  =  1.  (C.13) 

We  have  therefore  proven  that  our  eigen  functions  Me  indeed  orthonormal  for  both 
the  gravity-inertia  modes  and  the  geostrophic  modes. 


Appendix  D 

ON  THE  USE  OF  TEMPERTON’S  CFFT  ROUTINES 

In  this  appendix  the  complex  fast  Fourier  Transform  (CFFT)  routines  which  were 
used  in  the  Cartesian  coordinate  model  are  discussed. 

The  CCFT  routines  used  in  the  model  were  developed  by  Clive  Temperton  ( 1983a, b) 
and  can  be  obtained  through  NCARs  public  domain  software  library.  The  CFFT  package 
is  called  cfft99f.f  and  was  written  perform  a  number  of  simultaneous  complex  periodic 
Fourier  transforms.  The  code  was  written  to  run  efficiently  on  a  Cray  computer  through 
use  of  parallel  processing  techniques.  It  is,  however,  portable  to  other  machines. 

Although  FFT  and  CFFT  routines  are  well  documented  in  the  literature,  the  storage 
of  the  spectral  coefficients  is  explained  here  to  zud  in  understanding  the  model  code.  The 
cfft99f.f  routine  takes  a  complex  array  of  data  of  length  7  (0  to  7  —  1,  where  7  has  no 
prime  factor  greater  than  five)  and  returns  a  complex  array  of  equal  length  containing  the 
spectral  coefficients  of  original  data  (the  inverse  can  also  be  done).  What  is  important  to 
note  is  the  order  with  which  the  spectral  coefficients  are  returned.  The  first  zero  through 
7/ 2  storage  locations  in  the  array  contain  the  zero  through  7/2  spectral  coefficients.  The 
last  half  of  the  array  (7/2  +  1  to  7  —  1)  contains  the  i  —  I  spectral  coefficients.  That  is, 
C-k  —  C/_fc,  where  C  represents  any  spectral  coefficient.  As  mentioned  in  chapters  three 
and  four,  we  cannot  use  all  returned  wavenumber  if  we  wish  to  eliminate  aliasing  error.  If 
we  choose  to  keep  —  K  to  K  wavenumbers,  the  middle  portion  of  the  returned  array  must 
be  neglected. 

In  the  present  model,  we  must  perform  a  CFFT  in  two  dimensions,  both  x  and  y. 
This  results  in  a  two  dimensional  array  of  data.  The  storage  of  the  returned  spectral 
coefficients  follows  those  explained  above  in  both  dimensions.  In  two  dimensions,  we  must 
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Figure  D.l:  The  storage  of  the  spectral  coefficients  in  the  Cartesian  coordinate  model 
code.  The  shaded  area  represents  data  which  is  neglected  to  prevent  aliasing  error. 

neglect  a  cross  shaped  pattern  in  the  returned  array  (see  Fig.  D.l)  to  prevent  aliasing 
error. 

To  make  the  returned  arrays  more  convenient  for  usage,  we  rewrite  the  needed  data 
to  a  two  dimensional  array  of  length  (-K:K,-L:L).  When  we  perform  the  inverse  transform, 
we  then  simply  zero  out  the  wavenumbers  greater  than  |  K  |  and  |  L  |,  rewrite  the  array  in 
a  form  identical  to  the  original  array,  and  call  cfft99f.f  to  perform  the  inverse  transform. 

It  is  important  to  note  that  although  the  original  data  field  is  real,  we  write  the  data 
to  a  complex  array  to  perform  the  transform.  This  is  done  s^  cfft99f.f  can  be  used  in  both 
x  and  y,  thus  ensuring  compatibility  by  avoiding  the  use  of  a  real  FFT  routine  in  one 
direction  and  a  CFFT  in  the  other. 


END 


